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Abstract In this chapter of the e-book "Self-Organized Criticality Systems" we 
summarize some theoretical approaches to self-organized criticality (SOC) phenom- 
ena that involve percolation as an essential key ingredient. Scaling arguments, ran- 
dom walk models, linear-response theory, and fractional kinetic equations of the 
diffusion and relaxation type are presented on an equal footing with theoretical 
approaches of greater sophistication, such as the formalism of discrete Anderson 
nonlinear Schrodinger equation, Hamiltonian pseudochaos, conformal maps, and 
fractional derivative equations of the nonlinear Schrodinger and Ginzburg-Landau 
type. Several physical consequences are described which are relevant to transport 
processes in complex systems. It is shown that a state of self-organized criticality 
may be unstable against a bursting ("fishbone") mode when certain conditions are 
met. Finally we discuss SOC-associated phenomena, such as: self-organized turbu- 
lence in the Earth's magnetotail (in terms of the "Sakura" model), phase transitions 
in SOC systems, mixed SOC-coherent behavior, and periodic and auto-oscillatory 
patterns of behavior Applications of the above pertain to phenomena of magneto- 
spheric substorm, market crashes, and the global climate change and are also dis- 
cussed in some detail. Finally we address the frontiers in the field in association 
with the emerging projects in fusion research and space exploration. 
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1 The percolation problem 

The standard theory of percolation fV\ began with an attempt to make statistical 
predictions about the possibility for a fluid to filter through a random medium, pre- 
dictions that could be applied to a variety of physics problems, such as epidemic 
processes with and without immunization, the underground spread of pollution, and 
electrical discharges in thunderstorms. The phenomenon is characterized by a finite 
threshold, to be associated with a critical concentration of fractures, pores, or other 
sort of conducting channels in the medium, below which the spread is limited to a 
finite domain of ambient space, and is unlimited otherwise. The percolation prob- 
lem is relevant for a number of transport problems with threshold behavior as for 
instance Anderson localization [2| and hopping conduction in amorphous solids |3j. 
The percolation transition is perhaps the simplest phase transition-like phenomenon, 
with the macroscopic connectedness thought as a spontaneously occurring property, 
and the concentration of conducting elements as the control parameter (similar to 
the thermodynamical temperature) [4J. 

Site and bond percolation 

Given a periodic lattice, embedded in a li-dimensional Euclidean space, one can 
choose between two alternative formulations of the percolation problem: site and 
bond. The differences between site and bond percolation are actually very subtle and 
are manifest in a typically lower threshold for the bond problem. There also exists a 
hybrid, site-bond percolation due to Heermann and Stauffer |5|. In site percolation 
one assumes that the lattice sites are occupied at random with the probability p (and 
hence with the probability 1 — p are empty). A connected cluster is defined as a 
collection of all occupied sites that can communicate via the nearest-neighbor rule. 
In bond percolation, one thinks of clusters of connected conducting bonds instead. 
In this formulation all sites are initially occupied and bonds are occupied randomly 
with the probability p. Statistically, the p value decides on how big the connected 
clusters could be for the given topology of the lattice. The key point of the theory 
is the existence of a critical value, the percolation threshold pc (0 < pc < 1), above 
which the connected clusters span the entire lattice with the probability 1 . The crit- 
ical probability being smaller than 1 implies that the infinite clusters do not fill the 
ambient space yet. For p < pc, the percolation dies away exponentially. The thresh- 
old value is non-universal: it depends on the type of the percolation problem (site, 
bond, or hybrid); details of the lattice (cubic, diamond, triangle, etc.); as well as the 
ambient dimensionality d>l. The typical realizations of site and bond clusters on a 
square lattice are illustrated in Fig. 1. It is worth remarking that any point belonging 
to the infinite percolation cluster can be connected to the infinitely remote point via 
a connected escape path which lies everywhere on the cluster. For comprehensive 
reviews on percolation see, e.g., Refs. Ill |6l |7] HJ . 
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Fig. 1 Site vs bond percolation. Top: Site percolation problem, with circles representing occupied 
sites on a lattice. Left: A lattice system below the percolation threshold. Right: The same lattice 
system at the threshold of percolation, with a noticeably denser concentration of the occupied 
sites. Bottom: A similar picture for bond percolation. Connected (conducting) bonds are shown in 
blue color Left: An insulating state of the lattice below the percolation threshold. Right: Random 
distribution of conducting bonds at the threshold of conducting dc electricity. 

Percolation critical exponents ji, v, and jX 

Likewise to traditional critical phenomena, characterized by a scale-free statistics of 
the spontaneously occurring quantities, the geometry of connected clusters in vicin- 
ity of the percolation threshold is self-similar (fractal) 16] [51. As p — >■ pc, the perco- 
lation correlation (i.e., pair connectedness) length diverges as o= \p~ Pc\^^ ■ For 
p> the probability to belong to the infinite cluster is Poo (;?) °<: (p — Pc)^ 
whereas the dc conductivity behaves as Cdc °^ {p~ Pc)^ °^ t,^^!"^ . The critical ex- 
ponents j3, V, and /i are universal in that they do not depend on the type of the 
percolation problem, nor on details of the lattice. They do depend on the ambient 
dimensionality d, however, and their numerical values are known in all > 1. The 
"mass" of a connected cluster scales with its size as M(^) ^''p^^p) oc |Srf-^/v^ 
leading to a nontrivial Hausdorff dimension df = d — j5 /v. The latter expression is 
sometimes said to be "hyperuniversal" as it holds in all > 1. 
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Random walks on percolating clusters 

The problem of diffusion on fractals |[TOl [TTl [T2l has stirred considerable attention 
in the literature, especially, in terms of the random walk approach. If the random 
walker (an unbiased "ant") is put on a connected cluster at percolation, then the 
distance it travels after t time steps behaves as ll8l [T2]| 

(r2(f))c<f2/(2+e)^ 

The exponent Q is given by Q = {ji — j5)/v. Remark that the dependence here is no 
longer proportional to the time f, by contrast with uniform spaces. Thus, diffusion 
is anomalous. The exponent 9 describes topological characteristics of the fractal 
(such as connectivity, etc.) As such, it shows interesting invariance properties under 
smooth (diffeomorphic) maps of fractals [13J . It has two common names in the 
literature: the connectivity exponent and the index of anomalous diffusion. This 
ambiguity merely reflects that the diffusion is anomalous because fractals possess 
anomalous connectedness features as voids are present at all scales. When — >^ 0, 
normal (Fickian) diffusion is reinstalled. 

In a basic theory of percolation |6 1 it is shown that jU > )3 for connected clusters, 
implying that > 0. One sees that the mean-square displacement in Eq. ([T]l grows 
slower-than-linear with time. This slowing down of the transport occurs as a result 
of multiple trappings and delays of the diffusing particles in cycles, bottlenecks, and 
deadends of the fractal object on which the random motions concentrate. Note that 
the scaling law above holds as a single-cluster rule (the "ant" cannot jump between 
the clusters). Averaging over all clusters at percolation replaces Eq. ([T]) by 

{r'-{t)) oc f(2-/3/v)/(2+e) (2) 

for f <C Equation ^ has implications for the ac conductivity at "anomalous" 

frequency scales, O) 3> fQj- which the charge carriers move only on the 

fractal |12il14J. The various aspects of anomalous diffusion in fractal systems are 
summarized in the reviews ll8l [T3l[T5]| . 

The spectral fractal dimension and the Alexander-Orbach conjecture 

A hybrid parameter ds =2df/{2 + 9) is often referred to as the spectral, or fracton, 
dimension. It is so called because it represents the density of states for vibrational 
excitations in fractal networks termed fractons lfT6l[T7l . It also appears in the prob- 
ability of the random walker to return to the origin (oc f-'^»/2) The key differ- 
ence between the Hausdorff and spectral fractal dimensions lies in the fact that df 
is a purely structural characteristic of the fractal, whereas ds mirrors the dynami- 
cal properties (via the connectedness issues) such as wave excitation, diffusion, etc. 
Note that the spectral fractal dimension is not larger than its Hausdorff counterpart, 
provided that the connectivity exponent > 0. The value of ds can conveniently 
be considered as an effective fractional number of the degrees of freedom in fractal 



Percolation Models of Self-Organized Critical Phenomena 



5 



geometry, owing to the specific way it appears in the diffusion fT2l [TSl and wave- 
propagation problems ifTSl [TtI [T9| . 

In the past years there has been much excitement about the Alexander-Orbach 
(AO) conjecture |17| that the spectral fractal dimension is exactly 4/3 for percola- 
tion clusters in any ambient dimension d greater than 1. This conjecture is important 
as it relates the structural characteristics of the fractal, contained in dj, to the dy- 
namical characteristics, contained in Q. For d > 6, the AO conjecture was proven 
by Coniglio 1201 as a percolation problem on a Cayley tree (Bethe lattice). A Cay- 
ley tree is a graph without loops where each node contains the same number of 
branches (called the coordination number) 121]. In many ways, owing to its hierar- 
chical structure, a Cayley tree behaves as an infinite-dimensional space (its volume 
grows exponentially fast with the spatial scale). Not surprisingly, the percolation 
problem on a Cayley tree is regarded as a suitable model for mean-field percolation. 
For d < 6, the mean-field approach is invalidated as loops become important at all 
scales, thus impeding reduction to the trees. A great deal of effort has been invested 
to prove or disprove the AO conjecture in the lower embedding dimensions d < 6. 
It is now clear that in these dimensions the AO conjecture is not exact, nor does it 
generaUze to all statistical fractals as for instance to the backbones of percolation 
clusters l|8ll22l. Even so, the "true" values found for the spectral fractal dimension 
at criticality continue to be numerically surprisingly close to the original AO result 
ds =4-/3 for all d/ > 2 thus sustaining the conjecture. An interested reader may refer 
to the reviews L2.^. 

Percolation problem on the Riemann sphere 

It is both interesting and instructive to demonstrate how the spectral fractal dimen- 
sion may be obtained for threshold percolation on a plane {d ~ 2). The main idea 
here ||231 is to extend the plane by adding the point at infinity to it, then consider 
a stereographic projection of the infinite percolation cluster on the familiar from 
complex analysis Riemann sphere (see Fig. 2). Thus, the percolation problem will 
be compactified lfT3]| . as the point at infinity is mapped to the north pole. 

Observe that the percolation cluster spanning the plane implies that its stere- 
ographic image covers part of the surface of the sphere, the north pole included. 
Without loss of generality, we may assume that the south pole at which the sphere 
touches the plane belongs to the cluster. When considered on the Riemann sphere 
a percolating escape path to infinity originating from the south pole will be a sim- 
ple arc connecting the two poles, south to north. The key step is to notice that a 
solid angle at the base of this arc has a lower bound as posed by connectedness. 
Indeed this angle cannot be lesser than the angle at the base of half meridian. The 
latter is immediately seen to be equal to tt = \2n. Clearly, the percolation cluster 
itself is based on a solid angle not lesser than this. It is convenient to think of the 
number ds of the degrees of freedom as corresponding to an orthogonal basis of ds 
vectors ll23l . which span a fractional solid angle Q.^^ = dsTi'^^^^ /r{ds/2+ 1). Here, 
r denotes Euler's gamma function. Partial cases of this expression are, Q2 — 2?r 
for ds = 2 and Qt, — An for ds = 3. Thus, we expect that, for connected clusters. 
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Fig. 2 Stereographic projection of a percolating cluster on the Riemann sphere. The north pole, 
N, represents the point at infinity. When considered on the Riemann sphere a percolating escape 
path to infinity (blue line, with points A and B on it) originating from the south pole at which the 
sphere touches the plane is a simple arc connecting the two poles, south to north. A solid angle at 
the base of this arc has a lower bound, 7t, as dictated by connectedness. Then the mapping being 
conformal implies that at the threshold of percolation Qj^ = 7t, leading to ds = 1.327 ±0.001. 



^ ^> from which a lower bound on ds may be deduced by defining Q^^ — n. We 
associate this lower bound with the threshold of macroscopic connectedness (thresh- 
old of percolation). In equating X2j, to n we used that the stereographic projection 
being a conformal map is angle and circle preserving. Putting all the various pieces 
together, we have ifTSl l23l 

'^'r{dj2 + \)=^r{dj2)=''- 

Numerical solution shows that d^ = 1.327 ±0.001, remarkably close to, although 
slightly smaller than, 4/3. Rigorously speaking, this result disproves the AO conjec- 
ture in d = 2. Despite being this subtle, the observed deviation from 4/3 is impor- 
tant as it helps avoid the secular terms problem when applying a renormalization- 
group technique near the percolation point |7 8]. Note that the solution to Eq. Q 
has a remarkable meaning. It defines fractional dimensionality of a ball-like space 
seen from its center under the solid angle n. It is the dimensionality of this space, 
ds ~ 1.327, which permits percolation in terms of a connected escape path to in- 
finity. One sees that the percolation problem is essentially a topological problem. It 
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decides on macroscopic connectedness of random systems in terms of the number 
of the coupled degrees of freedom. As such, the percolation problem has impor- 
tant implications for the dynamics of complex systems, and self-organized critical 
systems as particular case, as it will be demonstrated shortly. 

Summary 

Summarizing, the percolation problem presents a non-trivial problem with scale- 
free behavior. It is related with phase transition-like phenomena as well as the fun- 
damental topology (via the connectedness issues). The percolation clusters provide 
a particularly clear example of statistical fractals in the limit ^ — > oo. The percolation 
indices j3, V, and ji are known in all ambient dimensions d > I and do not depend 
on details of the lattice nor the type of the percolation problem (universality). In ad- 
dition to basic science, the percolation problem is of practical importance as it offers 
a platform for the description of transport properties of disordered (random) media. 
In particular, diffusion and electrical conduction problems on percolation clusters 
have been widely studied and discussed in the literature (Refs. |I7][8][13]; references 
therein). 



2 The SOC hypothesis 

The challenge to understand fractals f24| and the 1// "noise" led Bak, Tang, and 
Wiesenfeld ||251 to introduce the concept of self-organized criticality, or SOC. The 
claim was that irreversible dynamics of systems with many coupled degrees of free- 
dom ("complex" systems) would naturally generate self-organization into a critical 
state without fine tuning of any external or control parameter(s). By analogy with 
traditional critical phenomena it was argued that in vicinity of the critical state there 
is universal behavior, robust with respect to variations of parameters and with re- 
spect to randomness, and that the system can be characterized by power-laws and a 
set of critical exponents. An impressive list of publication^have been produced in 
the attempt to prove or disprove the SOC hypothesis for the various systems. The 
phenomenon was demonstrated on a number of automated lattice models, or "sand- 
piles," displaying avalanche dynamics and scale invariance ll25l l26l l27l |28J . The 
various aspects of self-organized criticality dynamics have been reviewed by Bak 
ll29l . Jensen ll30l . Turcotte ||3T1 . Charbonneau et al. Il32l and Aschwanden ll33l . 

To qualify as SOC, the system must be open, be coupled with the exterior, and 
involve many interacting degrees of freedom. In addition, its dynamics must be 
thresholded and nonlinear, and the driving, or energy injection, rate must be very 
slow (infinitesimal). An important advance of SOC is the realization that fractals 
appear naturally through a self-organization process and that the corresponding crit- 
ical state is an attractor for the dynamics. In many ways the notion of SOC can 

' According to ISFs Web of Science the number of papers citing the seminal work in Ref. L25J is 
above three thousand. 
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be thought of as belonging to the nascent "science of complex systems" which ad- 
dresses the commonalities between apparently dissimilar natural, technological, and 
socio-economic phenomena, e.g., market crashes and climate disruptions |34|. De- 
spite its promising performance the SOC hypothesis is a subject of strong debate in 
the literature, and many issues related to it remain controversial or in demand for 
further investigation. 

SOC vs percolation 

Before we proceed with the main topics of this chapter, we would like to ad- 
dress the SOC hypothesis against the percolation problem discussed above. Indeed 
SOC shares with percolation the implications of threshold behavior and spatial self- 
similarity. An essential difference is that percolation is a purely geometrical model, 
whereas SOC involves, in addition, the temporal counterpart of the fractal, the 1 // 
noise ll35i . In many ways SOC is a spatio-temporal phenomenon where both spatial 
and temporal self-similarities are coupled and long-ranged. 

Another important aspect is that in percolation and other traditional critical phe- 
nomena, control parameters must be fine tuned to obtain criticality (thus the name 
"control"). In SOC phenomena, control parameters make part of the dynamical sys- 
tem instead: their values are defined dynamically as the system self-adjusts to ac- 
commodate the changing exterior conditions. It is in this sense that a SOC system 
is said to "automatically" (without a fine tuning of parameters) reach the critical 
state 0More so, we remark that there exists a "self-organized" formulation of some 
standard percolation processes such as the spread of deceases or forest fires. It was 
argued that their dynamics could be formulated so that they mimic SOC phenomena 
||36| . The characteristic feature here is that singularities at pc emerge not in distribu- 
tion of order parameters but of control parameters, making these phenomena look 
like SOC. This formulation is advantageous, as it leads to efficient numerical algo- 
rithms, allowing for a precise determination of the critical behavior, as for instance 
in models of self-organized critical directed percolation, with time interpreted as the 
preferred dimension 137,1 . 

The "guiding" mechanisms 

An important issue concerns the mechanisms that "guide" a system to criticality. 
These mechanisms are of two types. One type is associated with the application of 
an extremal principle that the dynamics should obey in order to satisfy the micro- 
scopic equations of motion. Examples of this type are the invasion percolation |38 1 
and Bak-Sneppen models ll39ll40l . In invasion percolatiorj^— introduced in physics 
by Wilkinson and Willemsen — the dynamics proceed along a path of least resis- 
tance under the action of capillary forces. Under the condition that the flow rate is 
infinitesimal the system finds its critical points that are stable and self-organized 

^ Often one says that a SOC system possesses no tunable control parameters, but that's all about 
the wording. 

' To be distinguished from ordinary percolation discussed above. 
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ll38l . The second type is associated with the operation of a nonhnear feedback be- 
tween system's dynamical parameters BTI as for instance a feedback of the order 
parameter on the control parameter(s) as discussed by Sornette B2ll43l . The Bak, 
Tang, and Wiesenfeld's (BTW) sandpile is a prominent example of this type. In 
sandpiles the unstable sand slides off to decrease the slope and reinstall stability, 
thus providing a feedback of the particle loss process on the dynamical state of the 
pile. In many ways nonlinearity is an essential key element to SOC phenomena as 
it ensures a steady state where the system is marginally stable against a disturbance 
MM- 



3 Going with the random walks: DPRW model 

The percolation problem when account is taken for a dynamical feedback mecha- 
nism offers a suitable platform to build toy-models of self-organized critical phe- 
nomena. Early attempts in this direction refer to the "dilution-by-hungry-ants" and 
the "thermal-fuse" models (with and without a healing) |42|. In what follows, 
we discuss a model p4l l45l . dubbed dynamic polarization random walk (DPRW) 
model, which combines the implication of a feedback mechanism with the idea of 
random walks on a fractal cluster at percolation. The model is formulated as a trans- 
port problem for electrically charged particles of different kinds|^ The advent of 
random walks in place of automated lattice redistribution rules makes it possible to 
calculate the frequency-dependent complex susceptibility of the dynamical system 
at SOC along with the memory (response) function and in the end to obtain the SOC 
critical exponents in terms of three percolation critical indices j5, jj., and v. This ap- 
proach paves the way for an analytical theory of SOC starting from the microscopic 
dynamical properties. One by-product of the random-walk model is a demonstra- 
tion 1 44, 45 1 that the relaxation of a supercritical system to SOC is of Mittag-Leffler 
type |46| (similarly to the relaxation in glassy systems and polymers). Indeed the 
Mittag-Leffler relaxation implies that behavior is multi-scale with a broad distri- 
bution of durations of relaxation events consistently with a description in terms of 
the fractional relaxation equation |47| and at odds with the well-known, Debye-like 
relaxation dynamics. 

Description of the model 

We consider a hypercubic J-dimensional {d>Y) lattice confined between two op- 
posite {d— l)-dimensional hyperplanes, which form a parallel-plate "capacitor" as 
shown in Fig. 3. The plate on the right-hand-side is earthed. Free charges are built 
by external forces on the capacitor's left plate. When a unit free charge is added to 
the capacitor the lattice responds by burning a unit "polarization" charge, which is 
an occupied site added at random to the lattice. When a unit free charge is removed 



* This electrical context of the model is non-crucial and can be relaxed 1451 . 
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from the capacitor a randomly chosen occupied site is converted into a "hole" site 
(missing occupied site). A hole will be deleted from the system (converted into 
empty site) if/when the corresponding free charge has reached (or been moved to) 
the ground. There is a limit, Qmax-> on the amount of the free charges the capacitor 
can store, and this is defined as 2max — epcN, where e is the elementary charge 
(e = —1), pc is the percolation threshold, and is the total number of sites across 
the lattice. Thus, the ability of the capacitor to store electric charges is limited to 
the occurrence of the infinite cluster at the percolation point. If, at any time, the 
above limit is exceeded, a double amounj^ of the free charges in excess of Q^ax 
will be removed from the capacitor and will be distributed between the sites of the 
infinite cluster with equal probability. The implication is that the capacitor leaks 
electric charges above the percolation point. This property reflects the onset of the 
dc conduction at the threshold of percolation. 



Sub-critical: P < P^. Super-critical: P > P^. 




Fig. 3 Dynamic polarization random walk (DPRW) model. Grey, blue, and azure particles show 
respectively the free charges, polarization charges, and holes. Left: System below the percolation 
point. Free charges are built by external forces on the capacitor's left plate. Right: System slightly 
above the percolation point, with an illustration of hopping activities on the lattice. Adapted from 
Ref. 145J. 

When a hole appears on the infinite cluster it causes an activation event with 
the following consequence: One of the nearest-neighbor occupied sites, which is a 
random choice, will deliver its charge content to the hole. The hole which has just 
received the polarization charge becomes ordinary occupied site, while the donor 
site becomes a hole. The newborn hole, in its turn, will cause further activation 
event at the location where it has occurred thus triggering off a chain reaction of 
redistribution of polarization charges. The chain reaction continues until the hole 
reaches the earthed plate where it is absorbed (converted into empty site). When 
a hole appears on a finite cluster it causes a chain reaction of activation events in 
much a similar way as on the infinite cluster but with one modification regarding 
the ending of the activation: The chain reaction stops if (i) likewise to the infinite 



This mimics non-zero inductance in the conduction process. 
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cluster case the hole reaches the earthed plate where it is converted into empty site, 
or if (ii) there are no more activities going on on the infinite cluster. In the latter 
case the finite cluster freezes in a "glassy" state with the still holes in it until either 
a new hole appears on the infinite cluster or one or more occupied site are added to 
the lattice by external forces. 

Random-walk hopping process 

Essentially, the holes interchange their position with the nearest-neighbor occupied 
sites, and it appears reasonable to model this process as interchange hopping process 
ll48l . In what follows we assume, following Ref. ifTOl . that there is a characteristic 
microscopic hopping time, which is taken to be unity, but more general hopping 
models can be obtained by introducing a distribution of waiting times between con- 
secutive steps of the hopping motion (continuous time random walks, or CTRW's) 
||49l l50 |. With the above assumption that the site acting as donor is a random choice 
the transport model is defined as random-walk hopping model. Similarly to the hole 
case, the free charges are assumed to behave as unbiased random walkers after their 
re-injection on the infinite cluster. They will hop at a constant rate between the 
nearest-neighbor occupied sites in random direction on the cluster on which they 
are initially placed until they reach the earthed plate where they get sink in the 
ground circuit (see Fig. 3). The holes act as the conducting sites for the motion of 
the free charges. The charged plate acts as a perfectly reflecting boundary. Hops to 
empty sites are forbidden. The latter condition limits the random walks to fractal 
geometry of the threshold percolation. 

Dynamical geometry of threshold percolation 

Overall, one can see that the system responds by chain reactions of random-walk 
hopping processes when it becomes slightly supercritical and it is quiescent other- 
wise. Excess free charges dissipating at the earthed plate provide a feedback mech- 
anism by which the system returns to the percolation point. There will be a slowly 
(as compared to hopping motions) evolving dynamical geometry of the threshold 
percolation resulting from the competition between the adding of occupied sites to 
the lattice and the charge-releasing chain reactions. Based on the quantitative anal- 
ysis below, we identify this state as a SOC state. This general picture based on the 
idea of a dynamic polarization response with random-walk hopping of the charge 
carriers might be called dynamic polarization random walk (DPRW) model. 

In the DPRW SOC model, chain reactions of the hopping motion acquire the role 
of "avalanches" in the traditional sandpiles. In the present analysis, we are interested 
in obtaining the critical exponents of the DPRW model by means of analytical the- 
ory. Numerical simulation of the DPRW dynamics is under way for comparison 
with the analytical predictions. By the time this chapter is being written, the char- 
acteristic signatures of multi-scale conductivity response of the dynamical system 
at criticality have been confirmed in the computer simulation model. In Fig. 4, we 
illustrate the existence of relaxation events of various sizes due to hole hopping on 
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a 10x10 square lattice with random distribution of the conducting nodes and the 
probability of site occupancy such as to mimic the percolation threshold and the 
conjectured SOC activities. 



DPRW model, 10x10 lattice 
1.5X10*1 




1000 2000 3000 4000 

time 



Fig. 4 Sizes of chain reactions due to hole hopping on a 10 x 10 square lattice. Each spike corre- 
sponds to a chain reaction of hole hopping on the system-scale conducting cluster, with the height 
proportional to the number of hops to absorption at the earthed plate. This simple numerical realiza- 
tion illustrates the existence of relaxation events of various sizes consistently with the implication 
of SOC. Adapted from Ref. gS). 



4 Linear-response theory 

Dynamics and orderings 

Starting from an empty lattice (no potential difference between the plates), by ran- 
domly adding occupied sites to it, one builds the fractal geometry of the random, or 
uncorrelated, percolation, characterized by three percolation critical exponents j3, 
V, and pL (connected clusters have fractal dimensionality df = d — P/v) |l4||6]|7l. 
Remark that the infinite percolation cluster, in the true sense of the wording, ex- 
ists only in the thermodynamic limit when the lattice itself is infinite. This limit 
arises because of the need to model the system-sized conducting clusters in terms of 
fractal geometry. In the absence of holes this percolation geometry is static (polar- 
ization charges can only move by exchanging their position with a hole) but when 
the holes appear on the lattice they cause local rearrangements in the distribution of 
the conducting sites. As a consequence, the conducting clusters on which the trans- 
port processes concentrate change their shape and their position in the real space. 
In the analysis of this section we shall require that the average number density of 
the holes be very small compared to the average number density of the polariza- 
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tion charges. The imphcation is that the system remains near the percolation point 
despite the slow evolution of the conducting clusters. Note that the lattice rules are 
such as to preserve the properties of the random percolation. In fact, no correlations 
are introduced in the distribution of the conducting sites at any step of the lattice 
update. 



Frequency -dependent conductivity and diffusion coefficients 

Given an input electric driving field E(f,r) the polarization response of the system 
is defined through 

F{t,r)= r°°Xit-t')W,r)dt', (4) 



where the response function x{f ~ f') is identically zero for f < f' as required by 
causality. We should stress that nonlocal integration over the space variable is not 
needed here in view of the local (nearest-neighbor) character of the lattice interac- 
tions. In a model in which the assumption of locality is relaxed as for instance in 
models permitting particle's Levy flights the integration over the space variable is 
expected to produce a physically nontrivial effect. We do not consider such models 
here. A Fourier transformed x{t) defines the frequency-dependent complex suscep- 
tibility of the system, x(g)). In a basic theory of polarization response one also 
introduces the frequency -dependent complex ac conductivity, ffac(fij), which is re- 



lated to xi('}) by the Kramers-Kronig integral [Eq. ( lOl below]. The dependence of 
the ac conductivity on frequency, specialized to the random walks on percolation 
systems fSjlll], is given by the scaling relation (7ac(o) °^ co^ , where the power ex- 
ponent (0 < 7] < 1) is expressible in terms of the percolation indices j3, V, and jj, 
as rj = jj. / {2v + fl ~ P) . We should stress that the scaling <7ac (c) °^ (0^ incorporates 
conductivity responses from all clusters at percolation including those finite. In the 
DPRW model these implications are matched by the mechanism of the hole con- 
duction permitting the polarization current on both infinite and finite clusters. The 
general linear-response theory expression iBTI for the conductivity (Tac(ci)) in terms 
of the mean-square displacement from the origin (r^(f)) is 

9 

ne 



9 

0-ac(«) = r^D{C0), (5) 



where 

—D{(o) = lim 

n^ £->-o+ 



{icof j\-'""e-" {r^{t))dt 



(6) 



with a constant depending on the dimensionality of the lattice and n and e the 
density and charge of the carriers, respectively. The function D{(£)) has the sense 
of the frequency-dependent diffusion coefficient ll52l . In the zero-frequency limit, 
Eq. (j5]) reproduces the well-known Einstein relation between the static diffusion co- 
efficient on the infinite cluster, Doo, and the dc conductivity, Odc — limco^o Caclo)- 
Note that the dc conductivity occurs only through the infinite cluster (p > pc), as 
opposed to the ac conductivity response, which occurs through both finite and in- 
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finite clusters. In what follows, we require that the frequency CO be large compared 
to the characteristic evolution frequency in the distribution of the conducting sites. 
Denoting this last frequency by O)*, we have, for the present orderings, 0)^0*. 
In this parameter range, D{(o) °^ (0^ . Consistently with the above definitions, the 
inverse frequency, 1 /(»*, is ordered as the characteristic diffusion time on the infi- 
nite cluster, T* ~ ^^/D(a)*). Note that this time will depend on O)* in accordance 
with Eq. (jsjl. Hence, O)* ~ I/t« ~ D(a)*)/^^, where ^ °^ \p — Pc\^^ is the diverg- 
ing percolation pair connectedness length; p is the probability of site occupancy; 
and pc is the percolation threshold. We have, at the margins of self-similar behav- 
ior, o= (o)*)'', implying that OJ* oc \p- pc\^^^^^^^'> . Observe that O)* ^ for 
p — >■ Pc. Remembering that there is a microscopic hopping time, which is taken to 
be unity, we assess the Kubo number 15311541 in vicinity of the SOC state as 

e,^lM^c<|;,-;,,|-v(l+'))/(l-')). (7) 

One sees that 2* — > oo in the limit p — > pc- The Kubo number is a suitable dimen- 
sionless parameter which quantifies how the evolution processes in the lattice com- 
pare with the microscopic hopping motions. The divergency of the Kubo number at 
criticality implies that there is a time scale separation: fast hopping motions vs in- 
finitesimal evolution change. In terms of the Kubo number (2* — °°), the diffusion 
coefficient D{(a^) becomes 

oc (o^Ql, (8) 

where we have introduced y= 1 — T] . This scaling law appears in models of anoma- 
lous diffusion by low-frequency turbulence fSS* '561 l57l l58l l59l . Special cases of 
Eq. ([8]l include the well-known Bohm scaling [60J, characterized by 7= 1, as well 
as the anomalous so-called "percolation" scaling (/Ri 0.7), dating back to diffusion- 
advection models of Isichenko and co-workers ll4l l6Tll62ll . Alternatively, the diffu- 
sion coefficient on a time varying fractal distribution, Eq. ([8]), can be deduced from 
the general scahng form ll59l 

{v\t))=e{ti^Af{th*), (9) 



where / is a scaling function, which interpolates between the initial-time power-law 
and flat asymptotic (f +°°) behavior: f{°°) = const. The form in Eq. (j9| is sim- 
ilar to that considered by Gefen et al. |12| for anomalous diffusion on percolation 
clusters (in their model, °^ ^^^^), and earlier by Straley lITTl . 

Power-law power spectral density 

By applying the Kramers-Kronig relations \mx{(o) °= (Jac(w)/fi^ and 

Rez(a)) - V.P. / , gac(ft)0 (10) 

J co'{co'-co) 
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it is found that xi^) °^ (i>^^ with 7 = 1 — 77. A Fourier transformed Eq. Q reads 
P((B,r) = ;i;(a))E(a),r). One can see that the power spectral density, 5(a)), of the 
system response to a white-noise perturbation, E(a),r) = 1, will be proportional to 
|;if(a))p. The end result reads: 



where a = 2(1 — 7]) = 2y. The conclusion is that the power spectral density in 
the DPRW model is given by an inverse power-law distribution, with the a value 
depending on scaling properties of the ac conductivity response. 

Stretched-exponential relaxation and the distribution of relaxation times 

Next, we obtain the distribution of relaxation times self-consistently. For this, as- 
sume that the system is slightly supercritical, then consider a charge density per- 
turbation, 5p(f,r), caused by the presence of either free charges or holes on the 
conducting clusters. "Slightly supercritical" means that the dependence of the ac 
conductivity response on frequency can, with good accuracy, be taken in the power- 
law form (7ac(ft)) °= 0)'^^ discussed above. The implication is that at adding 5p(f,r) 
to the conducting system at percolation we neglect the departure of the systems ge- 
ometric properties from pure self-similarity. Without loss in generality, we assume 
that the perturbation 5p{t.r) is created instantaneously at time f = 0. That means 
that the function 5p (f , r) = for f < for all r. The perturbation 5p (f, r) generates 
an electric field inhomogeneity, 5E(f,r) in accordance with Maxwell's equation 
V • 5E(f ,r) = 4;r5p(r,r). Consistently with the above discussion, we adopt that for 
t>Q the decay of 5p (f , r) is due to the spreading of charge-carrying particles (elec- 
trons and/or holes) via the random walks on the underlying fractal distribution. The 
polarization response to 5E(f,r) is given by 



where, as usual, x{t~t') = for f < f'. The density of relaxation currents is defined 
as the time derivative of 5P(f,r), i.e.. 



(11) 




(12) 




(13) 



The continuity implies that 



dt 



5p(f,r)+V. 




Xit-t')5E{t',r)dt' ^0. 



(14) 



Taking V- under the integral sign, then eliminating 5E(f,r) by means of Maxwell's 
equation V • 5E(r,r) = 47i:5p(f ,r), we find, with the self-consistent charge density. 
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5p{t,Y)+ATt x{t-t')5p{t',r)dt']^0. (15) 



d 
dt 

In writing Eqs. ^T4\ and ( fTSj ) we have also assumed that f > 0. We now integrate in 
Eq. ([T5]l to find 

5p{t,r)+47t j^\[t-t')5p{t',r)dt' ^g{r). (16) 

Here, the function g{r) is an arbitrary function of the position vector r, which ap- 
pears in the derivation as the constant of integration over time. Under the conditions 
X{t~t')=Q for r < t' and 5p(f,r) = for f < for all r, Eq. ( 16 1 reduces to 

5p(f,r)+47r f x[t -t')5p{t' ,r)dt' ^ g{r). (17) 



If we allow t +0, we find that for 7 > the integral term on the left-hand-side 
goes to zero (as 0= t^); 

lim r xit~t')5p(t\r)dt' ^ lim [' x(t ~t')5p(t' ,r)dt' ^0, (18) 

from which it is clear that g{r) = lim,^+o5p(f,r). We consider this last condition 
as the initial condition for the relaxation problem. Essentially the same condition 
holds in the limit 7^0, provided that lim,^+o is taken first. A Fourier transformed 
Eq. ( 17 I reads 

5p{co,k)+47txi(o)5pi(o,k)^g{k)/co, (19) 



where k is position vector in reciprocal space, and g(k) is the Fourier image of g{r). 
Writing the susceptibility as x(g)) = ''a)^''/47r with a time constant it is found 
that 

5p(o,k) = Vr^^lk). (20) 

The quantity has the sense of lifetime of a perturbation with wavelength A . We 
expect that oc at criticality, where z is a scaling exponent. A derivation of this 
scaling relation will be given shortly. Separating the variables, we write 5p(a),k) = 
(p(co)g(k), with 

which we consider as the relaxation function in the frequency domain. On inversion 
to the time domain, Eq. (21 1 generates the Mittag-Leffler function, Ey[—{t/Tx)'^, 
which has series expansion II46I l47l 

Thus, <p{t) = Ey[—{t/Txy]. One sees that the relaxation to SOC of a supercritical 
state is described by the Mittag-Leffler function Ey[~{t / Tx)'^^ ™d not by a simple 



= l/(» + T, ''w'-^), (21) 
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exponential function as for standard relaxation. We note in passing that the Mittag- 
Leffler function is the natural generalization of the exponential function. The latter 
is included as a special case 7=1. For times f ^ T;^, the Mittag-Leffler function, 
Eq. (j22]i, can be approximated by a stretched-exponential the so-called Kohlrausch- 
Williams-Watts (KWW) relaxation function | 



£y[-(f/T,)''] ^exp[-(f/T,)Vr(l + 7)], (23) 

which is often found empirically in various amorphous materials as for instance in 
many polymers and glass-like materials near the glass transition temperature (for 
reviews see Refs. (65] and f66l, and references therein). The KWW relaxation func- 
tion can conveniently be considered [67lj as a weighted average of the ordinary ex- 
ponential functions, each corresponding to a single relaxation event in the system: 

exp[-(f/T;i)Vr(l + 7)]= r e-'l^'wy{At)dAt. (24) 

Jo 

The weighting function Wy{At) is given by Eqs. (51d) and (55) of Ref. f67l where 
one replaces the exponent a with 7, the time constant T with T;^, and the variable pL 
with Tx/ At. In our notations: 

Wy{At) = {Tx/Ae)Ly,^y{Tx/At), (25) 



where Ly._i is the Levy distribution function with skewness —1 (e.g., Ref. f68i). 
Assuming a long-wavelength perturbation (i.e., the parameter X being much longer 
than the microscopic lattice distance: A ^ 1 ), and setting ix/ Af:^ 1 , we can further 
approximate the Levy distribution Ly _i by the Pareto inverse-power distribution. 
This gives Ly^^\ [ix/At) oc {■Zx/At)-^^+y'^ leading to a pure power-law distribution 
of relaxation times, consistently with the wisdom of SOC: 

Wy{At)ocAt-^At^+'f^At-'^. (26) 

These distributions were earlier conjectured for SOC f2E\. Our conclusion so far 



is that the relaxations are multi-scale, in accordance with Eq. (24i, and their du- 
rations are power-law distributed. The distribution is heavy-tailed in the sense that 
/ dAtWy{At) oc — !• 00 for Xx °°- 

Consistency check 

In a basic theory of dielectric relaxation one writes the frequency-dependent com- 
plex dielectric parameter as Il67ll69l 

e{co)-loc- r^^e""dt, (27) 
Jo at 



where (p{t) is the relaxation function that describes the decay of polarization after 
the polarizing electric field has been stepped down or removed instantaneously. In 
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the DPRW model, a step-down type electric field occurs as a consequence of re- 
injection of the free charges to the infinite cluster. The ensuing relaxation dynamics 
are mimicked by the chain reactions of hole hopping which act as to properly redis- 
tribute the polarization charges across the lattice. Based on the above analysis, we 



identify the relaxation function in Eq. (27 i with the Mittag-Leffler function to yield 
(p{t) ~ Ey[—{t/Txy]. In vicinity of the critical state, because the upper limit on Tx 
diverges, we can, moreover, replace Ey[—{t/Txy] with exp[— (f/T;L)''/F(l + 7)] for 
(almost) all < f < 00. Thus, for p pc, (p{t) ~ exp[-{t /r^y /r{l +y)]. Integrat- 



ing by parts in Eq. ( 27 1, after a simple algebra one obtains 

£{co)-locl-sV{s)+isQ{s), (28) 

where s = coTx is dimensionless frequency, and Q{s) and V(.s) are the Levy definite 
integrals: 

Q(s)~ / exp(— M'')cos(Mi)iiM, (29) 
Jo 

Wis) — / exp (—u"^) sin {us)du. (30) 
Jo 

In the parameter range of multi-scale relaxation response, T^/^t ^ 1' (O'^A. ^ 1- the 
following series expansions of the Levy integrals hold |67|: 

QW-Lc-ir'^'^-f^- (31) 

,^1 s"r+^ r(« + l) 2 

VW=L(-l)"^ U"'^^ os^. (32) 
,^0 ^"^ r(n+l) 2 



From Eqs. (31 1 and (32 1 one can see that the expansion of Q{s) starts from a term 
which is proportional to and so does the expansion of V(.s) — 1 /s. Hence, up 

to higher order terms, £{(o) — 1 0= s^^. Given this, one applies the Kramers-Kronig 
relations sQ{s) °^ aac{(o)/(0 and 

1 _^,V(.V) - V.P. f^j^ rCJac(fi)') (33) 

to find the scaling of the ac conduction coefficient to be (7ac(o) °^ fi)'^''. By com- 
paring this with the above expression C7ac(o) °^ co^ one reiterates that 7 = 1 — 
consistently with the distribution of durations of relaxation events, Eq. ([26]). 



Fractional relaxation and diffusion equations 

As was shown by Glockle and Nonnenmacher llTOl . the Mittag-Leffler function. 



Eq. ( 22 1, is the solution of the fractional relaxation equation 

rI^ = -oD;->(0, (34) 
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where 

is a fractional time the so-called Riemann-Liouville derivative BtIITTII . Partial cases 
of this derivative are the unity operator for 7 ^ 1 and d/dt for 7 — >^ 0. It is no- 
ticed, following Ref. |72|, that the Mittag-Leffler function Ey[-'{t /T^y] describes 
the relaxation toward equilibrium of particles governed by the fractional diffusion 
equation, or FDE 

|:P(r,r)=oD^^v2p(f,r), (36) 

where P{t,r) is the probability density of finding a particle (random walker) at time 
t at point r, and the Laplacian operator stands for the local (nearest-neighbor) char- 
acter of the lattice interactions. 



Derivation of the fractional diffusion equation 



It is instructive to obtain the fractional diffusion equation, Eq. ( 36 1, directly from 



the DPRW relaxation model. For this, let us introduce the electrostatic potential, 
5^{t,r), corresponding to the electric field inhomogeneity, 5E(f ,r) = — V5^'(f,r). 



Upon substituted into Eq. ( 14 1 



X{t-t')V^8<P{t',r)dt'. 



(37) 



In the vicinity of self-organized critical state, we can represent the total charge den- 
sity, p (f , r), as a sum of "unperturbed" or background density, = \e\pc, and a per- 
turbation, 5p{t,r), describing the deviation from criticality: p{t,r) = Pc + 5p{t,r). 
To obtain the dependence of p{t,r), we use the effective-medium approximation 
II73I . a standard technique for calculating average physical properties of many -body 
systems. The idea is to think of the particles as embedded into an "effective" poten- 
tial, 50{t,r), where they will be Boltzmann-distributed in accordance with 



p(f,r)=p,exp[e50(f,r)/r]. 



(38) 



Self-consistently, one requires that, on the average, the embedding in the effective 
medium has the same overall property as the effective medium itself ||48] |73l. In 



writing Eq. (38 1 we took into account that the perturbation, 5p(f,r), is due to neg- 
atively charged particles (electrons and/or holes). The normalization condition is 
defined through p(f,r) Pc for 50{t,r) — > 0. In the above, the parameter T has 
the sense of "thermodynamic temperature," associated with the random motion of 
particles on a conducting cluster. For |e5^'(f,r)|/r <C 1, expanding the exponential 
function on the right of Eq. (38 1, one finds 5p(f,r) w {pce^ /T)5't>{t,r). Thus, for 



small perturbations (high temperatures), 5p{t,r) is proportional to 5^'(f,r), as it 



should. Eliminating 5<P{t,r) in Eq. (37 1, we have 
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^5p(f,r) = I l^\{t-t'){T/p,e^)V^5p{t\r)dt'. (39) 

The memory function, x{t),is obtained as Fourier inversion of xico) = '^(O^^ /An, 
yielding x{t)°<^t'^^'^ . Under the conditions x{t -t')=Qfoxt <t' and 5p (f , r) = for 



f < for all r, the improper integration in Eq. ( 39 1 can be performed in the limits 
from to t. Collecting all dimensional and numerical parameters in one effective 
"diffusion coefficient," Ay ^ {T /4-7tpce^)T^ ^, we can write 



which is an equivalent form of Eq. p6f . The fractional diffusion equation, Eq. ( 40 1, 
can be thought of as deriving from the generalized "Fick's law" 

^^(^'•■^-rk^i'lT^^'^^''^^''' '''' 

or 5j(f,r) = —qDI ^Ay,Vdp{t,r), which is readily deduced from Eq. ( [l3| in the 
effective-medium approximation. Equation ( |4T] i can equivalently be obtained |59| 
from the general scaling law for anomalous diffusion on percolation systems. A 
derivation using the scheme of continuous time random walks (CTRW's) can be 



found in Refs. 1,74. ,75 J . We should stress that the non-Markovian nature of Eqs. (40 1 



and (41 1, together with the fractional relaxation equation, Eq. ( 34 1, accounts for the 
long-time memory effects in SOC phenomena, where one believes li42l the time 
evolution exhibits long tails and infinite correlation scale. 

The occurrence of the fractional diffusion equation, Eq. ( [40| , might be inter- 
preted, with the aid of the proposed SOC model, in favor of considering SOC as one 
important case for fractional kinetics [76] . The concept of fractional kinetics enters 
different areas of research, such as turbulent transport in plasmas and fluids, parti- 
cle dynamics in potential fields, quantum optics, and many others. This subject is 
summarized in comprehensive reviews ll47ll77ll78J . In many ways equations built on 
fractional derivatives offer an elegant and powerful tool to describe anomalous trans- 
port in complex systems |77|. There is an insightful connection with a generalized 
master equation formalism along with a mathematically convenient way for calcu- 
lating transport moments as well as solving initial and boundary value problems 
BtI ITS I . The fundamental solution or Green's function of the fractional Eq. ( [36| is 
evidenced in Table 1 of Ref . [78J . 



Dispersion-relation exponent, Hurst exponent, and the X-exponent 

In sandpile SOC models, one is interested in how the lifetime of an activation cluster 
scales with its size |27|. In the DPRW model by activation cluster one means a 
connected cluster of activated sites. An occupied site is said "activated" if it has 
become a hole or if it contains a free charge. Clearly, activation clusters can only 
exist above the percolation threshold. Note that activation clusters are subsets of 
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the underlying conducting cluster of polarization charges. The notion of activation 
cluster is but a visualization of the charge density inhomogeneity 5p(f ,r) in terms 
of a connected distribution of activated sites. Activation clusters decay because the 
constituent charged particles (holes and/or free charges) diffuse away via the random 
walks. 

Consider an isotropic activation cluster composed of the free particles. (The na- 
ture of the particles does not matter here — the hole case is just similar) It is assumed 
for convenience, without loss of generality, that each site of the activation cluster 
contains only one particle. Thus, the number density of the free particles inside the 
activation area is equal to one. It steps down to zero just outside. If the microscopic 
lattice distance is a (a = 1), then there is a unit density gradient across the boundary 
of the activation cluster looking inside. Because of this gradient the activation clus- 
ter will be loosing particles on the average. A particle that has crossed the boundary 
against the direction of the gradient is considered lost from the cluster As the par- 
ticles dissipate, the location of the density pedestal shifts inward with speed u. The 
local flux density of those particles leaving per second the activation area is just the 
gradient times the local diffusion coefficient. The latter depends on frequency of 
the relaxation process as D{(o) °^ (0^ in accordance with Eq. (jsj). If / is the current 
size of the cluster, then the corresponding relaxation frequency is O) ~ m/Z. Using 
this, the frequency dependence of the diffusion coefficient can be translated into 
the corresponding /-dependence, the result being D{1) oc Balancing the rate of 
decay of the cluster with the outward flux of the particles we write dl/dt oc — 
Integrating this simple equation over time from t — to t = and over I from 
/ = A to / = one finds the dispersion relation °^ X' with z=l + Tj=2 — 7. The 
persistency |9 | of relaxation is measured by the Hurst exponent H, which is related 
to our z via H — l/z- Lastly, the T-exponent, which defines the distribution of the 
particle flows caused by a single chain reaction, or activation-cluster size distribu- 
tion II27I . is obtained from Eq. (5) of Ref. |26|, where one replaces g with a, and 
the fractal dimension D with the fractal dimension of the infinite percolation cluster, 
df = d — P /v. The end result is T = 3 — az/df. Note that the T values in Refs. Ii26l 
and differ by 1. 

Occurrence frequency energy distribution and the ^-exponent 

It is convenient to think of the activation clusters as containing a certain amount 
of "energy" which is released when the comprising free particles dissipate to the 
boundaries. Using here that the electric charge of the free particles is a conserved 
quantity, we may associate the energy content of the activation clusters with their 
electric charge content. The latter — accounting that each site of the activation clus- 
ter contains only one particle — will be proportional to the fractal volume ~ A^/, 
with A the size of the cluster. More so, the energy confinement time, Tx, can be 

zldf 

obtained as the activation cluster lifetime, enabling Xx °^ A' 0= £^ , where the dis- 
persion relation has been used. Differentiating in Eq. ( [26] l, one arrives at the inverse- 
power distribution 

N{ex)<^dwyiex)/dex<^e^\ (42) 
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where 6 = 1 + rjz/df is a power-law slope and df = d — j5/v in the "hyperuniversal" 
fractal dimension. We interpret this distribution as the familiar from applications 
Il32l l79l occurrence frequency energy distribution. Indeed the distribution in Eq. (42 1 
might find its significance in the statistics of solar flares and gamma bursts, provided 
that the underlying physics processes refer to SOC. We estimate the 6 values shortly. 



Values of the critical exponents 

Using known estimates ||6l|4]|71 of the percolation indices j3, V, and we could 
evaluate the critical exponents of the DPRW model in all ambient dimensions d>\. 
The results of this evaluation, summarized in Table 1, are in good agreement with 
the reported numerical values from the traditional sandpiles (for c/ = 2, z « 1 .29, 
T « 2.0; for li = 3, z w 1.7, T w 2.33) |26| and earlier theoretical predictions (for 
c/ = 2, z 4/3, T = 2; for 3, z = 5/3, T = 7/3) |27|. We consider this conformity 
as a manifestation of the universality class of the model. For d ~ the model 
reproduces the exponents of mean-field SOC lISOl . 

In many ways, the DPRW approach to SOC offers a simple yet relevant lattice 
model for dielectric relaxation phenomena in systems with spatial disorder One by- 
product of this approach is the case for stretched-exponential the KWW relaxation 
function, Eq. (23 i, which is often found empirically in various amorphous materials 
as for instance in many polymers and glass-like materials near the glass transition 
temperature (for reviews see Refs. |65 66|; references therein). In this connection, 
we observe that the model gives values of the exponent 7 (for d = 2, y^i 0.66; for 
d = 3, ysa 0.4) in good agreement with the typical experimental results (the / value 
between 0.3 and 0.8) |65 671 HT] 182]. This observation supports the hypothesis 
ll83l[84l that dielectrics exhibiting stretched exponential relaxations are in a state of 
self-organized criticality. 

More so, the DPRW model gives a Hurst exponent (for d — 2,H 0.75; for d — 
3, // « 0.6) consistently with the reported narrow range of variation of H as observed 
in different magnetic confinement systems (Hurst exponent varying between H w 
0.62 and 0.75) (85] |86l |87l |88l . In this connection, it is worth remarking that SOC 
behavior of the bulk plasma transport is expected to be a characteristic of higher- 
power plasma discharges in the so-called low confinement regime 1 85 1. 

With respect to the occurrence frequency energy distribution, Eq. (42 1, the model 
predicts that 6 = 1 in one dimension (for d =\,r\ =Q) and 6 = 3/2 in the mean- 
field limit (for t/>6, z = 2, Tj = l, and df =4). These results are exact. Also, one 
finds, approximately, 6 w 1.24 for li = 2 and 6 w 1.4 for d — 1> (Hausdorff fractal 
dimensions df = 91/48 and df w 2.5, accordingly). These theoretical predictions 
are in close agreement with the predictions of Aschwanden Ii89il who found values of 
6 = 1 .0 for c/ = 1 ; B = 1 .28 for d = 2; and 6 = 1.5 for d = 3 from a statistical fractal- 
diffusive avalanche model of a slowly-driven self-organized criticality system (Ch. 
2.2.2, this volume. In the notation of Ref |891, 6 = a^). The latter model assumes 
that the avalanche size grows as a diffusive random walk, implying that z = 2 in 
all dimensions d = 1,2,3. In our model, diffusive random walks are recovered in 
the mean-field limit only and in relatively high ambient dimensions that are not 
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lesser than 6, while the z exponent is taken to be non-diffusive in general. Even 
so, this does not seem to affect the 6 value significantly, such that our mean-field 
result, 13 = 3/2, almost precisely coincides with the result of Aschwanden |89| in 
three dimensions. Earlier, Litvinenko |90| has suggested that the distribution of flare 
energies is characterized by a power- law with the slope 6 = 3/2 independently of 
the ambient dimensionality d > \.Y\t modeled an avalanching process on a tree 
without loops, thus giving rise to this value. In this context, we should stress that 
the effect of loops can be abandoned only in the high dimensions d >6, permitting 
a mean-field description |7 8|. All in all, the exponent 6 = 3/2 agrees well with 
the reported slopes of the occurrence frequency energy distribution for solar flares 
(around —1.5 to —1.8) ll79l 1911 l92l . demonstrating that the observed power-law 
distribution of flare energy release is well reproduced under the assumption that the 
solar corona operates as a self-organized criticality system 1,32., 33., 89.. 93.. 94J . 

Table 1 Critical exponents of the DPRW model. Exponents appearing in statistical distributions 
as for instance inverse power-law power spectral density distribution are summarized in the lower 
part of the table. The mean-field results, holding for J > 6, are collected as J = oo. Input parameters 
are the percolation indices j8, v, and fi (4)|6j,7 |. 
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5 Discussion 

General 

Apart from details of the mathematical formalism, the DPRW model is actually 
quite simple. The main points are as follows. A lattice site can either be empty 
or occupied. An occupied site is interpreted as polarization charge. The equilib- 
rium concentration of the polarization charges depends on the potential difference 
between the plates. When the potential difference changes the lattice occupancy 
parameter adjusts. A dynamical mechanism for this uses holes. The holes are just 
missing polarization charges. They are important key elements to the model as they 
provide a mechanism for the polarization current in the system. Beside holes, the 
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free charges are introduced. The free charges, too, carry electric current whose very 
specific role in the model is just to control the potential difference between the 
plates. The changing amount of the free charges in the system has effect on the lat- 
tice occupancy parameter Nonlinearly, it affects the conductivity of the lattice. This 
nonlinear twist provides a dynamical feedback by which the system is stabilized 
at the state of critical percolation. In many ways the proposed model is but a sim- 
ple lattice model for dielectric relaxation in a self-adjusting disordered medium. It is 
perhaps the simplest model which accounts for the whole set of relaxation processes 
including the hole conduction. 

It is worth assessing the advantages and disadvantages of the DPRW approach to 
SOC. In terms of advantages, the electric nature of the model greatly facilitates the 
analytical theory: Not only does it permit to quantify the microscopic lattice rules in 
terms of the frequency-dependent complex ac conductivity, the use of the Kramers- 
Kronig relation in Eq. (lOi makes it possible to directly obtain the susceptibility 
function by integrating the conductivity response. As a result, the exponents z, 7, a, 
and H are expressible in terms of only one parameter, the exponent of ac conduction 
r\. The latter is obtained as a simple function of the percolation indices j3, V, and jU. 

With respect to disadvantages, the model is seemingly different from the tradi- 
tional approaches to SOC based on cellular automation (CA) and its integration in 
the existing family of SOC models might be a matter of debate. Even so, the idea 
of the random walks on a self-organized percolation system as a simplified yet rele- 
vant model for SOC constitutes a significant appeal: First, it relies on the established 
mathematical formalism of the random walks lIU [Ts) 07] [tU whose advance on the 
SOC problem is theoretically very beneficial. Second, it offers a clear connection 
to studies outside the conventional SOC paradigm as for instance to transport of 
mass and charge in disordered media BSl [STl l52l [95 1. Instead, the traditional CA 
type models are complicated by a poor analytical description of the microscopic 
transport mechanisms and their basic physics appreciation is at times uneasy. 



The role of random walks 

It is theoretically important to note that the dielectric context of the considered 
model, apart from offering a convenient platform for the analytical theory, is not 
however essential for the SOC phenomena. Indeed the DPRW model could be de- 
fined in terms of diffusion processes for neutral particles of different kinds. A formal 
reason for this is the equivalence |51 52| of the frequency-dependent electrical con- 
ductivity problem and the frequency-dependent diffusion problem, specific to hop- 
ping conduction. The crucial element to the model is the assumption of the random 
walks, not the nature of the particles. 

The possible generalizations of the DPRW model correspond to biased random 
walks of the free charges in the direction of the potential drop and/or inclusion 
of a second critical threshold p^c > Pc above which the random walk dynamics 
might change to a biased motion. We consider those generalizations obvious as they 
mainly intend to modify the value of the exponent T] in a certain parameter range, 
while the basic physical picture of SOC will remain essentially the same. 
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Universality class 

The final point to be addressed here concerns the issue of universality class. We 
take notice of the fact that the DPRW SOC model uses the charitable redistribution 
rule ll37ll to propagate the activities, likewise to the traditional BTW sandpile \25] 
or similar ll26i, 27 1. That means that an active site always loses its content to the 
neighbors. The charitable rule is to be distinguished from the neutral rule, when 
each of 2(i + 1 sites involved in redistribution gets an unbiased random share of the 
transported quantity. Models using the neutral rule often fall in the universality class 
of directed percolation (DP) and are characterized by appreciably larger values of 
the dynamic exponents (for d = 2,z~ 1.73 ±0.05) |37|. Based on this evidence, we 
suggest that the DPRW model belongs to the same universality class as the BTW 
sandpile, and not to the DP universality class, consistently with the values of the 
critical exponents collected in Table 1 . 



6 "Sakura" model 

The DPRW model of SOC can be extended so that it includes the phenomena of self- 
organized "turbulence" in Earth's geomagnetic tail, discussed in Refs. |96 97, 98|. 
The main idea here is that, when the processes of magnetic reconnection stretch the 
magnetotail beyond a certain limit, the cross-tail electric current system of Earth's 
dipole-like magnetic field is destabilized, since it crucially requires the presence of 
a regular component of the geomagnetic field normal to the current sheet plane. 
Then the magnetotail spontaneously evolves into a far from equilibrium dynamical 
"turbulent" state, where it responds to changes in the tail current intensity and the 
varying dawn-dusk potential difference in terms of turbulent perturbation electric 
currents and magnetic field fluctuations. It was argued that the transport of electric 
charge across the magnetotail was due to heavier plasma species, the ions, whose 
regular orbits (transient or Speiser, weakly trapped, trapped) were essentially de- 
stroyed in the absence of stabilizing normal magnetic field and that the steady state 
which is self-organized corresponded to a strongly shaped electric current system 
and to a highly inhomogeneous, multi-scale magnetic fluctuation pattern (Fig. 5). 

This model for the coupled turbulent perturbation electric currents and magnetic 
field fluctuations has come to be known as "Sakura" model, after its presentation 
(96 1 at the Chapman Conference held in Kanazawa, Japan, November 5-9, 1996. 
Applications of the Sakura model pertain to the phenomena of tail current disrup- 
tion in substorm regions of the near-Earth tail |98, 99l as well as to the explana- 
tion of permanent presence of magnetic field fluctuations in the distant magnetotail 
inilSTJ, as suggested by the GEOTAIL measurements ifTOOlfToTI . The important 
role of the ion component of the plasma in driving cross-tail current instability was 
addressed in Refs. II102II1031I104I . where one also finds an analysis of the satellite 
observational data. 
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Fig. 5 "Sakura" model. Top left: Schematic illustration for the Earth's dynamic magnetosphere 
interacting with its solar wind drive. Top right: The various orbit types for themial ions in the mag- 
netotail current sheet with finite regular component of the normal field. Adapted from Ref. 1 104|. 
When the processes of magnetic reconnection stretch the magnetotail beyond a certain limit, the 
particle regular orbital motion is essentially destabilized and the cross-tail electric current system 
spontaneously evolves into a far from equilibrium dynamical "turbulent" state where it responds 
to changes in the tail current intensity and the varying dawn-dusk potential difference in terms 
of turbulent perturbation electric currents and magnetic field fluctuations. Bottom left: Original 
model promotion as a blossoming sakura tree: Its "leaves" represent the magnetic field fluctua- 
tions; its "branches," the perturbation electric currents. Bottom right: Plan view of the turbulent 
current sheet. Cross-tail electric currents are organized in highly branched, very inhomogeneous 
conducting patterns with the topology of a fractal network at percolation (red color). The magnetic 
field fluctuations are shown as chunks of different sizes (violet color). They scatter the momentum 
of current-carrying particles, the ions, and are electromagnetically related with the perturbation 
electric current intensity by means of Maxwell's equation, V x 5B(?,r) = (4;r/c)5j(?,r). Note 
that the magnetic field fluctuations, 5B{f,r), can be thought as analog polarization charges in the 
DPRW model of SOC. 



Denoting the perturbation tail current density and magnetic fluctuation field by 
respectively 5j(f,r) and 5B(r,r), we write 

Vx5B(f,r) = — 5j(f,r), (43) 

c 

where, under the assumptions of locality and linear conductivity response. 
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5j(f,r) = j a{t~t')dE{t',r)dt'. (44) 

Here, 5E(f ,r) is the perturbation electric field in the current sheet plane. The mem- 
ory function, a(f ), is obtained as Fourier inversion of the frequency-dependent com- 
plex ac conductivity, a^dio), where the ac dependence is due to the turbulent nature 
of the conducting domain. It is understood that the time varying magnetic pertur- 
bation generates a time and spatially varying electric field because of Faraday's 
law. An important feature which arises in this induction process is gradual heating 
and energization of the plasma, leading to the occurrence of slowly decaying, high- 
energy non-thermal wings in the particle energy distribution function lfT3l 1 10511 1061 . 
Often in the magnetospheric plasma research those distributions with wings are 



modeled by non-thermal the so-called "kappa" distributions 11071 which interpolate 
between the initial low-energy exponential forms and the asymptotic inverse power- 
law behavior The significance of the "kappa" distributions lies in the fact 1 108 1 that 
they appear as canonical distributions in the non-extensive thermodynamics due to 
Tsallis (Ch. 2.3.5, this volume) 1 109|. In the present analysis we shall assume, how- 
ever, that all fluctuation frequencies are so slow that the effect of the inductive field 



can be neglected and we omit, consequently, the inductive term in Eq. ( 43 1 above 



To this end, performing a Fourier transform of Eq. (43 i in space and time, we have 



4-71 

kx 5B(o),k) = — 5j(a),k), (45) 

c 



where k is the wave vector of the perturbation. Simultaneously, from Eq. (44 1 we 



find 5j(a),k) = (7ac(fi^)5E(a),k). In vicinity of the current instability threshold, con- 



sidering the low-frequency limit of Eq. (43 i, we may think of the fluctuations as of 



collection of plane waves, characterized by a linear dispersion relation O) = k • u, 
where the phase velocity, u, does not depend on k. Given that the input perturbing 



electric field is uncorrected white noise: 5E(a),k) ~ 1, with the aid of Eq. (45 1 one 
sees that the power spectral density of the magnetic fluctuation field, |5B(a),k)p, 
will be proportional to S{(o) °= |<7ac(co) /cop. At this point, if one assumes, following 
Refs. 19611981 . that the dynamical "turbulent" state is characterized by fractal geom- 
etry of the threshold percolation, one may exploit the scaling relation C7ac(fi)) °^ co^ 
to obtain S{(o) °^ CO^", with a = 2(1 — ?]), consistently with the DPRW result. Uti- 
lizing the percolation estimate above (for d = 2,7] ^ 0.34), one finds that a « 1.3. 
This theoretical prediction compares well against the reported a values in the lower- 
frequency part of the magnetic fluctuation spectrum (below a turnover or knee fre- 
quency posed by the unstable tearing modes: see Fig. 6) BIOOI 11021 II 101. Thus, 
behavior is self-similar in the self-organization domain. Remark that the power 
spectral density of the magnetic fluctuation field, |5B(a),k)p, corresponds with the 



power spectral density of the dynamic polarization response, Eq. (111. We should 
stress that the Sakura model leads to a smaller spectral index (a sa 1.3) than Kol- 
mogorov's theoretical value for fluid turbulence (a = 5/3) as well as Kraichnan's 
theoretical value for magnetohydrodynamic plasma turbulence (a = 3/2). This last 
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observation addresses the significance of "self-organized" magnetic fluctuation tur- 
bulence as opposed to ordinary (fluid-like) turbulence. 



Fig. 6 The typical spectra of magnetic field fluctuations in the near-Earth stretched and thinned 
magnetotail prior to tail current disruption. Left: INTERBALL-1 observations. Courtesy of prof. 
Lev Zelenyi. Right: A spectrum observed by the Charge Composition Explorer of the Active Mag- 
netospheric Particle Tracer Explorers (AMPTE) satellite, with an emphasis on the August 28, 1986, 
current disruption event. Adapted from Ref. 1 102|. A distinctive feature of these spectra is the knee 
around a characteristic frequency posed by the unstable tearing modes. The spectrum is flatter be- 
low the knee frequency and is noticeably steeper just above it. The processes of self-organization 
to a critical state correspond to the flatter counterpart of the spectrum and to frequencies lesser 
than the knee frequency: in practice, lesser than ~ 5 • 10^^ Hz, although the exact bound, in real 
data, may not be that certain. Associated spectrum is a power-law with the typical slope ~ — 1.3 
(in log-log plot). In the limit of very low frequencies, the observed signals cross over to a white 
noise, as they should (this is due to finite system size effects and/or the natural limitations of the 
observational time series), so that the spectra are flat. The spectral properties of the fluctuations just 
above the knee (slope ~ —2.4) involve the processes of convection of magnetic turbulence struc- 
tures with decelerated solar wind velocity along the magnetotail. These processes were analyzed in 
Refs. 1 98 99 1 and are not considered here. All in all, the signatures of self-organization to a critical 
state are to be expected in the intermediate frequency range, comprised between those frequencies 
where the spectra are white noise-like and the knee frequency. An approximate position of the knee 
is marked by vertical dashed line. 

Our conclusion so far is that the power spectral density of the magnetic fluctua- 
tion field is given by an inverse power-law and that the behavior corresponds with 
the prediction of the DPRW SOC model. This result suggests that micro-processes 
of self-organization of electric currents and magnetic field fluctuations in the Earth's 
stretched and thinned magnetotail are governed by SOC llT3ll99l . Remark that the 
above assumption of fractality and self-similar behavior is validated through the 
direct analysis of time series of the satellite observed magnetic field fluctuations, 
obtained in situ in the substorm regions of the near-Earth tail and associated with 
the phenomena of tail current instability during the substorm growth phase II102I . 




Frequency [Hz] 
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7 Beyond the linear theories: DANSE formalism 

Theoretical approaches discussed so far were based on a linear-response theory and 
on fractional generalizations of the diffusion and relaxation equations, Eqs. ([34| 



and Eq. (40 1, which are linear by construction. Nonlinearities were contained in 
fractal geometry of the dynamical system at percolation and implicitly in the vari- 
ous "avalanche" exponents and the fractional indices of time differentiation. In the 
present analysis, this paradigm of "linear dynamics in a nonlinear medium" will 
be relaxed. Rather, a more general theoretical picture will be drawn, in which the 
dynamical and structural nonlinearities are twisted, and for which one might pro- 
pose the formula "nonlinear dynamics in a nonlinear medium." Very specifically, 
we intend to demonstrate that the phenomena of SOC — at least those belonging 
to the universality class of the DPRW model — may be cast in the mathematical 
formalism of discrete Anderson nonlinear Schrodinger equation (DANSE), which 
invokes a random potential for lattice interactions, and in which the strength of non- 
linearity, being an inherent part of the model description, is determined dynamically 
as the system self-adjusts and evolves to criticality in response to external forcing. 
Most previous theoretical studies of SOC have neglected the possibility of describ- 
ing the lattice interactions in terms of a random potential field, and have focused, 
consequently, on the microscopic redistribution rules for the dynamics. We believe 
that this approach is unnecessarily restrictive and has left out the important physics 
results. 



The roadmap 

As is akeady mentioned above, in the DPRW model the critical state is made self- 
organized via the mechanisms of hole-hopping by which the system responds to 
the fluctuating potential difference on the capacitor. We should stress that the holes 
redistribute the polarization charges in a way as to preserve the properties of the ran- 
dom percolation. They change the shape and the folding of the percolation clusters 
in the ambient configuration space, but not the random character in the distribu- 
tion of the conducting sites. A confirmation of these ideas can be obtained if one 
considers holes as "excitations" of self-organized critical state mediating the lattice 
activities. With this interpretation in mind, the transport problem for holes can be 
formulated as a transport problem for the hole wave function in a random potential 
field^The latter problem, in its turn, can be studied as part of the general problem 
of transport of waves in disordered media, with that very specific element that the 
medium is nonlinear and its properties are coupled with the wave process itself — 
which, too, can be nonlinear 

We consider the problem of dynamical localization of the hole wave function in 
the framework of nonlinear Schrodinger equation (NLSE) with a random potential 
term on a lattice. The NLSE was derived for a variety of physical systems under 
some approximations [ 1 1 1 J . Recently, it has been rigorously established that, for a 



* Essentially the same approach applies to the electrons. 
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large variety of physical conditions and details of interaction, the NLSE (also known 
as the Gross-Pi taevsky equation) is exact in the thermodynamic limit 111 121 . An 
important feature which arises in this approach is competition between randomness 
and nonlinearity. As we shall see, this competition has a significant effect on SOC 
problem. Under most general conditions, we can expect that, when the nonlinearity 
is sufficiently small, the random properties dominate, giving rise to the phenomena 
of Anderson localization of the excitations. That means that diffusion is suppressed 
and a wave packet that is initially localized will not spread to infinity |2|. 

So, what happens with the increasing strength of nonlinearity? The question is far 
than trivial, after the simulation-motivated suggestion 1 11 3 , 11 4 1 that a weak nonlin- 
earity can destroy localization above some level, giving rise to unlimited spreading 
of the wave field along the lattice, despite the existing disorder. The dynamics of the 
spreading has remained a matter of debate 1 1 14, 115 1 16 11 ITlfl 181 . 

In a recent investigation of NLSE with disorder, it has been theoretically found 
111 191 that destruction of Anderson localization in the presence of nonlinearity is a 
critical phenomenon — similar to a percolation transition in random lattices. De- 
localization occurs spontaneously when the strength of nonlinearity goes above a 
certain limit. Below that limit, the field is localized similarly to the linear case. In 
the analysis of this section, we bring these ideas in contact with the physics of SOC 
and we consider a situation in which the strength of nonlinearity is determined dy- 
namically in terms of time depending probability of site occupancy as the system 
fluctuates near the critical point. It is this very specific nonlinear twist with the dy- 
namical state of the lattice, which captures the essential key signatures due to SOC, 
and which as is suggested below enables to cast the SOC problem into the mathe- 
matical formalism of a discrete NLSE with disorder. The effect this twist has on the 
dynamics is that the localization-delocalization transition is turned self-organized, 
occurring exactly at the state of critical percolation. 

With the aid of our DPRW lattice model, we can essentially simplify the analysis 
if we consider that the transport in vicinity of the critical state occurs as a result 
of hole hopping between the nearest-neighbor sites, occupied by the polarization 
charges. As the percolation threshold is approached, we can envisage clusters of the 
polarization charges as the effective (random) medium, acting as a random poten- 
tial on the hole wave function. Following Anderson |2 |, we adopt the tight binding 
description for the hopping processes. Consequently, we introduce a Hamiltonian, 
paving the way to consider the transport problem for SOC as essentially a Hamil- 
tonian problem. This approach poses a theoretical challenge, as it aims to connect 
SOC with first-principle models. 

DANSE equation 

Focusing on the hopping motions on a lattice, we consider a variant of discrete 
Anderson nonlinear Schrodinger equation (DANSE) with randomness 

ih^ = HlWu + Clr«l V«, (46) 
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where 

HlWu = e„Wn+V{-^„+i + v/„_i), (47) 

Hl is the Hamiltonian of the Hnear problem in the tight-binding approximation 
i/A„ is the hole wave function; HiWn describes hopping-like transitions between the 
nearest-neighbor sites on a lattice; and ^|i//„pv/„ accounts for the generic nonlin- 
earity of the wave process. In the above, C, characterizes the strength of nonlin- 
earity]^ on-site energies e„ are randomly distributed with zero mean across a finite 
energy range; V is hopping matrix element; and the total probability is normal- 
ized to L„ I V'nP — 1- More general models can be obtained by replacing |v^„P with 
1 1^,1 1', for arbitrary r > 0. We do not consider such models here. In what follows, 
h = I for simplicity. When ^ 0, all eigenstates are exponentially localized 0. 
In the absence of randomness, DANSE (46i is completely integrable. Adhering to 
the effective-medium description, we aim to comprehend the spreading of the hole 
wave function under the action of nonlinear term in the limit t — >■ +0°. 



Coupled nonlinear oscillators 

It is useful to expand the hole wave function using an orthogonal basis of the eigen- 
states, 0„„,, of the Anderson Hamiltonian, Hl, yielding 

r„ =l^a,„(f)0„.m (m=l,2,...). (48) 

m 

"Orthogonal" means that Y.n^nm'Pn-k = ^m,k, where 5„,jt is Kronecker's delta and 
star denotes complex conjugate. The total probability being equal to 1 implies that 
UnVnV^n = L;« (0 (0 = 1- For the nonlinear equation, Eq. (46 1, the depen- 
dence of the expansion coefficients, C7,„(f), is found to b^ 

idk- 0}k(yk = C, ^ V<;_mi,m2,'n3^mi^m2^'«3' ^^^^ 

where dot denotes time differentiation; cOk {k= 1,2,...) are the eigenfrequencies of 
Hl; and the amplitudes Vit.mi.nn.ms are given by 

Vk ^n,mi- (50) 



In deriving Eq. (49 1 we took into account that HL<^n,k = (Ok^l^n.k- Equations (49 1 cor- 
respond to a system of coupled nonlinear oscillators with the Hamiltonian 

H = Ho+Hi„u Ha = Y,(^kO*kOk, (51) 

k 



^ We assume that the nonlinearity is repulsive > 0), implying that it favors the spreading of 
the wave field. In the opposite case of attractive nonlinearity < 0), solitons are typically found 

fTBllTTTl. 

^ Hint: substitute Eq. |48| into DANSE |46| , then multiply the both sides by and sum over n, 
remembering that the modes are orthogonal. 
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^int — 2 ^ ^k,m[,m2,'n3^k^lnl(^m2^l"3^ (52) 

k,mi ,m2,m3 

Thus we have translated the hopping problem for the hole wave function into the 
interaction problem for coupled nonlinear oscillators on a lattice. In the above, Hq 
is the Hamiltonian of non-interacting harmonic oscillators and Hi„t is the interaction 
Hamiltonianl^Each nonlinear oscillator with the Hamiltonian 

c 

hk = cokolok + -ijVkMXk'^k'^kOk'^k (53) 
and the equation of motion 

iOk - (Ok(yk - CVk.k..k.k<^kCTkCJk = (54) 

represents one nonlinear eigenstate in the system — identified by its wave number k, 
unperturbed frequency (Ok, and nonlinear frequency shift AcOk — C^k.k.k.k<^k<^k ■ Non- 
diagonal elements Vk^„ti.m2,m3 characterize couplings between each four eigenstates 
with wave numbers k, mi, m2, and mj. It is understood that the excitation of each 
eigenstate is not other than the spreading of the wave field in wave number space. 
Resonances occur between the eigenfrequencies (0/^ and the frequencies posed by 
the nonlinear interaction terms. We havq^ 

0>k = O/Hi - «m2 + <»m3 • (55) 

When the resonances happen to overlap, a phase trajectory may occasionally switch 
from one resonance to another As Chirikov realized II120L any overlap of reso- 
nances will introduce a random element to the dynamics along with some transport 
in phase space. Applying this argument to DANSE ( [46] l, one sees that destruction 
of Anderson localization is limited to a set of resonances in a Hamiltonian system 



of coupled nonlinear oscillators, Eqs. (51 1 and (52i, permitting a connected escape 
path to infinity. 



Chaotic vs pseudochaotic dynamics 

At this point, the focus is on topology of the random motions in phase space. We 
address an idealized situation first where the overlapping resonances densely fill the 
space. This is the familiar fully developed chaos, a regime that has been widely 
studied and discussed in the literature (e.g., Refs. 1 121 122 1). A concern raised over 
this regime when applied to Eqs. (jSTJ and ( [52| ) comes from the fact that it requires a 
diverging free energy reservoir in systems with a large number of interacting degrees 
of freedom such as SOC systems. Yet, developed chaos offers a simple toy-model 
for the transport as it corresponds with a well-understood, diffusive behavior. 



' We include self-interactions into H\„f . 

Conditions for nonlinear resonance are obtained by accounting for the nonlinear frequency shift. 
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A more general, as well as more intricate, situation occurs when the random mo- 
tions coexist along with regular (KAM regime) dynamics. If one takes this idea to 
its extreme, one ends up with the general problem of transport along separatrices of 
dynamical systems 1 123 1. This problem constitutes a fascinating nonlinear problem 
that has as much appeal to the mathematician as to the physicist. An original im- 
portant promotion of this problem to "large" (spatially extended) systems is due to 
Chirikov and Vecheslavov II124I . 

This type of problem occurs for slow frequencies. One finds ll58l |59l that 
resonance-overlap conditions are satisfied along the "percolating" orbits or sepa- 
ratrices of the random potential where the orbital periods diverge. The available 
phase space for the random dynamics can be very "narrow" in that case. In large 
systems, the set of separatrices can moreover be geometrically very complex and 
strongly shaped. Often it can be envisaged as a fractal network at percolation as for 
instance in random fields with sign-symmetry Il3ll4l ll25l . 

There is a fundamental difference between the above two transport regimes 
(chaotic vs near-separatrix). The former regime is associated with an exponential 
loss of correlation permitting a Fokker-Planck description in the limit f The 
latter regime when considered for large systems is associated with an algebraic loss 
of correlation instead 117711126111271 . implying that the correlation time is infinite. 
There is no a conventional Fokker-Planck equation here, unless generalized to frac- 
tional derivatives |74, 75 1, nor the familiar Markovian property (i.e., that the dy- 
namics are memoryless). On the contrary, there is an interesting interplay between 
randomness, fractality, and congelation, which is manifest in the fact that all Lya- 
punov exponents vanish in the thermodynamic limit, despite that the dynamics are 
intrinsically random |59|. 

This situation of random non-chaotic dynamics with zero Lyapunov exponents, 
being in fact very general, has come to be known as "pseudochaos" ||771|128II129L 
One may think of pseudochaos as occurring "at the edge" of stochasticity and chaos, 
thus separating fully developed chaos from domains with regular motions. There is 
a growing belief that the concept of pseudochaos offers the natural mathematical 
platform to obtain the fractional kinetic equations from first principles 1771 . 

In section 3 above it was argued that the phenomena of SOC could be described 
by fractional kinetics, which is a suitable and powerful formalism for long-range 
correlated behavior Here, we lay more stress on this argument by proposing that 
SOC processes are as a matter of fact pseudochaotic. In this way of thinking one 
naturally bridges the concepts of fractional kinetics, non-Markovian transport, and 
SOC. Then the inherent "edge" character of pseudochaotic dynamics can be related 
with what one believes is the threshold nature of SOC processes Il25l . Support to 
this suggestion can be found in the results of Refs. [i44l l45l [59l . 

Nearest-neighbor rule 

This idea of "edge" behavior brings us to a model 01 1911 where each nonlinear oscil- 
lator, Eq. ( [53] l, can only communicate with the rest of the wave field via a nearest- 
neighbor rule. Indeed this is the marginal regime yet permitting an escape path to 
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infinity. We associate this regime with the onset of delocalization. Clearly, the num- 
ber of coupling links is minimized in that case. Note that the nearest-neighbor rule 
guarantees that all interactions are local, this being an essential key element to SOC. 
When summing on the right-hand-side, the only combinations to be kept are, for the 
reasons of symmetry, akG^Git and aA._iC7j'(7t+i. We have 

i&k ~ (OkOk = i;,Vkakalak + 2!;v^ak-ialak+i, (56) 

where we have also denoted for simpHcity Vk — Vk,k,k,k and Vj^ — Vk,k-\.k.k+i - Equa- 



tions ( [56) l define an infinite {k — 1,2,...) chain of coupled nonlinear oscillators 
where all couplings are local (nearest-neighbor-like). The interaction Hamiltonian 



in Eq. (52 1 is simplified to 



r 

^ k k 



Pseudochaotic dynamics on a Cayley tree 

We are now in position to introduce a simple lattice model for the transport. The key 



step is to observe that Eqs. (56 1 can be mapped on a Cayley tree where each node 
is connected to c = 3 neighbors (here, c is the coordination number). The mapping 
is defined as follows. A node with coordinate k represents a nonlinear eigenstate, 
or nonlinear oscillator with the equation of motion ( [54) l. There are exactly c = 3 
branches at each node: one that we consider ingoing represents the complex ampli- 
tude <7^, and the other two, the outgoing branches, represent the complex amplitudes 
Ok-\ and Ok+\ respectively. These settings are schematically illustrated in Fig. 7. 

A Cayley tree being by its definition | |2T| a hierarchical graph offers a suitable 
geometric model for infinite-dimensional spaces. We think of this graph as embed- 
ded into phase space of the Hamiltonian system of coupled nonlinear oscillators. 



Eqs. (51 1 and (52 1. In the thermodynamic limit ^max — > in place of a Cayley tree, 
one uses the notion of a Bethe latticep] Setting /cniax — ^ we suppose that each 
node of the Bethe lattice hosts a nonlinear oscillator, Eq. ( [54| . The bonds of the 
lattice, in their turn, can conduct oscillatory processes to their neighbors as a result 
of the interactions present. 

Next, we assume that each oscillator can be in a chaotic ("dephased") state with 
the probability p (and hence, in a regular state with the probability 1 — p). The 
p value being smaller than 1 implies that the domains of random motions occupy 
only a fraction of the lattice nodes. Whether an oscillator is dephased is decided 
by Chirikov's resonance-overlap condition — which may or may not be matched 
on node k. We believe II 1241 that in systems with many coupled degrees of freedom 
each such "decision" is essentially a matter of the probability. The choice is random. 
Focusing on the p value, we consider system-average nonlinear frequency shift 



" A Bethe lattice is an infinite version of the Cayley tree. To this end, a purist might prefer to say 
"bond" in place of "branch," but that's all about the terminology. 
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Fig. 7 Mapping Eqs. p6| on a Cayley tree. Each node represents a nonlinear eigenstate, or nonlin- 
ear oscillator with the equation of motion itTj. — (Oj-Oj- — j8Vj.j(.j(- j-Oi.trJ'cFj. = 0. Blue nodes represent 
oscillators in a chaotic ("dephased") state. Black nodes represent oscillators in regular state. One 
ingoing and two outgoing branches on node k(k= 1,2,...) represent respectively the complex am- 
plitudes (Tj*, Oi— 1, and Gk+i- Structures that are not explicitly shown are beyond the dashed lines. 
Adapted from Ref fTT9l . 



zia)NL = C(lr«l')4„ (58) 

as an effective "temperature" of nonlinear interaction. It is this "temperature" that 
rules over the excitation of the various resonant "levels" in the system. With this 
interpretation in mind, one writes p as the Boltzmann factor 

p = exp(-5c(j/ziC(JNL), (59) 

where 5a) is the characteristic energy gap between the resonances. Expanding 
over the basis of linearly localized modes, it is found that 

{Wn?-)An = ^Y. L C„„0«.m2<,O-'"2- (60) 

n m\.m2 

The summation here is performed with the use of orthogonality of the basis modes. 
Combining with Eq. (58 1, 

r 

4tONL = -^rCT>m- (61) 
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The sum over m is easily seen to be equal to 1 due to the conservation of the prob- 
abiUty. Thus, /IOnl ~ C,/An. When the field is spread over An states, the distance 
between the resonant frequencies behaves as 5(x) ^ I /An. We normalize units in 



Eq. (46 1 to have 5(0 — 1/ An exactly. One sees that 

p = exp(-l/C). (62) 

For the vanishing ^ ^ 0, the Boltzmann factor p — > 0, implying that all oscillators 
are in regular state. In the opposite regime of ^ — > oo, p — > 1 . That means that all 
oscillators are dephased and that the random motions span the entire lattice. 

There is a critical concentration, pc, of dephased oscillators permitting an escape 
path to infinity for the first time. This critical concentration is not other than the 
percolation threshold on a Cayley tree. In a basic theory of percolation it is found 
that Pc = l/(c — 1) (e.g., Ref. ETIl '). This is an exact result. For c = 3, = 1 /2. 
We associate the critical value pc = 1/2 with the onset of transport in the DANSE 



model, Eq. (46 1. When translated into the ^ values the threshold condition reads 

C. = l/ln(c-l). (63) 

Setting c = 3, we have — l/ln2 w 1.4427. This value defines the critical strength 
of nonlinearity that destroys the Anderson localization. For the ^ values smaller 
than this, the localization persists, despite that the problem is nonlinear. When ^ > 
l/ln2, the localization is lost and the wave field spreads to infinity. 

Our conclusion so far is that destruction of Anderson localization is a thresholded 
phenomenon, which can be described as a percolation transition in a system of de- 
phased oscillators on a Cayley tree (Bethe lattice). Delocalization occurs when the 
strength of nonlinearity, mathematically related with the concentration of dephased 
oscillators, Eq. ( [62) l, exceeds a certain critical level. The critical point is exactly at 
Cc = l/ln2. 



Making delocalization transition self-organized 



With the recognition that, according to Eq. (62i, the nonlinearity parameter, ^, is 



twisted with the probability of site occupancy, p, the formalism of discrete Ander- 



son nonlinear Schrodinger equation, Eq. (46 1, allows for a representation, which 
makes it possible to include the phenomena of SOC. The main idea here is to think 
of ^ as a fluctuating parameter, which is defined dynamically as the system evolves 
to percolation, and whose value is decided "on-the-fly" by the actual state of the 
lattice. Then the ^ value need not fine tuned to its critical value, i^^, in order for the 
delocalization to occur, but it rather emerges as an attracting (singular) point as the 
original system of interacting charge-particles self-adjusts to percolation. It is this 
feature which leads to the dynamical rule of advancing the hole wave function and to 
a "self-organized" formulation of the localization-delocalization transition. We reit- 
erate that the critical strength of nonlinearity which destroys Anderson localization 
is expressible in terms of the percolation threshold as = — l/lnpc- This last result 
suggests that behavior be non-perturbative in vicinity of the critical point. Theoret- 
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ically, this observation is very important as it elucidates the nature of thresholded 
phenomena including SOC. It is also manifest in the pseudochaotic character of the 
lattice activities, implying that the correlations persist despite that the microscopic 
dynamics are inherently random 1,59. .77. .119. ,129J . 



Asymptotic spreading of the hole wave function 

Let us now obtain second moments for the threshold spreading of the wave-field. 
This task is essentially simplified if one visualizes the transport of the hole wave 
function as a random walk over a system of dephased oscillators. For p — > pc, this 
system is self-similar, i.e., fractal. It is this fractal distribution of dephased oscil- 
lators which, according to Bak, Tang, and Wiesenfeld |25|, conducts "the noise 
signal. . . through infinite distances" just above the marginally stable state. We con- 
sider this distribution as a percolation cluster on a Cayley tree. The fractal geometry 
of this cluster is fully characterized by the Hausdorff dimension df —4 and the in- 
dex of anomalous diffusion, 0—4 (e.g., Refs. |7, 8|). Note that the spectral fractal 
dimension is exactly 4/3 in this limit, consistenly with the original AO result ifTTll . 
From Eq. ([T|l one obtains 

((Zin)2(f)) f^+oo. (64) 

This behavior is asymptotic in the thermodynamic limit. Remark that the Hausdorff 



dimension being equal to 4 matches with the implication of Eqs. ( 49 1 and ( 50 1 where 
the coefficients Vjr.mi./n^.His supposed to run over 4-dimensional subsets of the 
ambient mapping space. Indeed it is the overlap integral of four Anderson eigen- 



modes, Eq. (50 1, that decides on dimensionality of subsets of phase space where 
the transport processes concentrate. When the nearest-neighbor rule is applied, this 
overlap structure is singled out for the dynamics. Under the condition that the struc- 
ture is critical, i.e., "at the edge" of permitting a path to infinity, the support for the 
transport is reduced to a percolation cluster on a Bethe lattice — characterized, along 
with the above value of the Hausdorff dimension, by the very specific connectivity 
exponent, 0—4. The end result is2/(2 + 0) = l/3. 



Summary 

We have shown that the Anderson localization in disordered media can be lost in the 
presence of a weak nonlinearity and that the phenomenon is critical (thresholded). 
That means that there is a critical strength of nonlinearity above which the wave field 
turns to unlimited spreading. Below that limit, the field is localized similarly to the 
linear case. We have discussed this localization-delocalization transition as a perco- 
lation transition on the separatrix system of discrete nonlinear Schrodinger equation 
with disorder. This problem is solved exactly on a Bethe lattice. A threshold for de- 
localization is found to be = l/ln2 w 1.4427. For the ^ values smaller than this, 
the localization persists, despite that the problem is nonlinear. Support for this type 
of behavior can be found in the results of Ref. IJBOJ . More so, a "self-organized" 
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formulation of the localization-delocalization transition has been obtained on the ba- 
sis of DANSE ( |46) l by defining the nonlinearity parameter on-the-fly, thus offering 
a fertile playground to describe the phenomena of SOC. In vicinity of the delocal- 
ization point the spreading of the wave field is subdiffusive, with second moments 
that grow with time as a powerlaw oc f V3 for t +00. This regime bears signatures 
enabling to associate it with the onset of "weak" transport fTsTl of Alfven eigen- 
modes (AEs) in the vicinity of marginal stability of magnetic confinement systems. 
In this respect, we note that the characteristic aspects of sandpile physics involving 
SOC have, in the AE transport case, been discussed in Refs. 1 1321 1 1331 . 



8 The two faces of nonlinearity: Instability of SOC 

Often when SOC systems are said to be "nonlinear" one refers to the operation of a 
feedback mechanism |42, 41] ensuring that the control parameters need not be fine 
tuned explicitly to obtain criticality. It is this feedback which stabilizes the system 
at the state of marginal stability, or the SOC state, as opposed to traditional critical 
phenomena, which do require a tuning. The comprehension of nonlinear feedback 
mechanism leads to another face of self-organization, the existence of a bursting in- 
stability of SOC 1 44, 45 1, which occurs in the parameter range of excessive external 
forcing, and for which we suggest the name "fishbone-like" instability, by analogy 
with some bursting (internal-kink) instabilities in magnetic confinement systems 
II134I . The implication is that nonlinearity can play either stabilizing as in the ideal 
SOC regime or destabilizing as in the regime of overdriving, role depending on the 
strength of interaction with the exterior. This section is aimed to discuss this topic 
in more detail with the aid of the DPRW model. 

Instability cycle 

As the probabiUty of site occupancy p approaches the percolation threshold pc, 
the pair connectedness length diverges as ^ 0= \p — Pc\^^- For p > p^, the longest 
relaxation time in the system is °^ {p ~ Pc) '^ ™d the dynamic susceptibility 
behaves as x °^ {p ~ Pc)~^^'^ ■ In the DPRW model p as a function of time is deter- 
mined dynamically by the competing charge deposition and loss processes. That is, 
dp/dt = Z+ — Z_, where Z+ is the net deposition rate of the free charges on the 
capacitor's left plate, and Z is the particle loss rate. The net deposition rate, or the 
driving rate, is the control parameter of the model: It takes a given value. The particle 
loss rate is obtained as electric current in the ground circuit, i.e., Z = 19{p — Pmm), 
where pmin is the lower limit of variation of p. Note that Z is due to the free 
particles leaving the system through the earthed plate. The Heaviside Q function 
indicates that the lattice can release charges only if/when p > Pmm- We expect that 
Pmin lies close to, although somewhat lower than, the percolation threshold pc- This 
is because the conducting cluster can still loose its charge content to the ground cir- 
cuit even in the absence of a connecting path to the charged plate. The dynamics of 
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/ can be estimated from dl/dt — ±1 /T^, where the upper sign corresponds to the 
relaxation process in the lattice. Flitting all the various pieces together, we write 

dp/dt=Z+-W{p-pmm), (65) 
dl/dt = WI\p - Pc^ sign{p - /?,), (66) 

where sign(/? — Pc) — +1 (—1) for p > pc (p < pc) is the sign-function, and W is 
a numerical coefficient. Equations ( [65] l and ( [66| define a simple system of equa- 
tions for two cross-talking variables, the lattice occupancy per site, and the parti- 
cle loss current. These model equations are perhaps the simplest nonlinear equa- 
tions describing the generic fueling-storage-release cycle in driven, dissipative, 
thresholded dynamical systems. An examination of these equations shows that the 
dynamics are periodic (auto-oscillatory), with the peak value of electric current 
/max ^ W^(Pmax - Pcf^^^ ■ Here, pmax (Pmax > Pc) is the Upper limit of the p vari- 
ation. Note that ly^ax ^ for p^^,^ pc as expected. The auto-oscillatory motions 
signify that the pure SOC state is destabilized and that the systems phase trajecto- 
ries enter the supercritical parameter range. When /?niax ^ 1, the periodic dynamics 
acquire a sharp, bursting character. The bursts half-duration (or half-width) equals 
A ~ (1/W)(;?max ^ Pc)^^^ ■ Eliminating the distance to the critical state one obtains 
the scaling relation /^ax °^ 4 where b ~ [zV + l)/zv. The period between the 
bursts is found to be ©f, ~ {pc — Pmm)/Z+. A pure SOC state with no superimposed 
periodic bursts arises when ©h — °°. This implies that Z+ — > for p — > pc- Thus, 
criticality requires the vanishing of Z+, in agreement with the result of Ref. lISOll . 
The critical state is stable when @t,^T^. We have 

{Pc-P^m)/Z+ » - \p-Pct'''. (67) 
This is satisfied when Z+ — > faster than 

Z+max b-^'c h'^- (68) 



This limit exceeded, the system turns to auto-oscillate around the percolation point 
with a period dictated by the net deposition rate of the polarization charges. The 
physics origin of this auto-oscillatory motion lies in the fact that the changing 
amount of the free particles provides a feedback on the lattice occupancy parameter. 
It is due to this feedback relation that the DPRW system operates as a self-adjusting, 
intrinsically nonlinear dynamical system. Whether or not this feedback will excite 
the instability depends on how the characteristic driving time compares to the char- 
acteristic relaxation time. Indeed, focus on the stability condition in Eq. (67 1. The 
system being stable at the percolation point requires that the relaxation time due to 
the random walks be short compared to the characteristic driving time, 1/Z+. 
In this parameter range, any occasional charge density perturbations will dissipate 
via the random walks before input conducting sites are again introduced. When the 
percolation point is approached, because the time scale oc \p — pc\ diverges, it 
is essential that the system be driven infinitesimally slowly to remain at a pure SOC 
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State. Instability occurs when the relaxation processes operate on a longer time scale 
than the driving processes. In this regime, the system accumulates the polarization 
chargesj^ whereas to remain at criticality it would get rid of them. The accumu- 
lation of the polarization charges has direct effect on the conductivity between the 
plates, which steps up with the lattice overshooting the percolation threshold. When 
Pmax — > 1, the system can be thought of as facing the typical conditions of elec- 
trostatic discharge in the regime of short-circuit. It should be emphasized that the 
feedback mechanism does a two-fold job: (i) it stabilizes the system at the state 
of critical percolation in a regime when the driving rate is infinitesimal; and (ii) it 
excites a cross-talk between the conductivity and the lattice occupancy parameters 
when the driving rate is faster than the relaxation rate. 

In the parameter range in which the strength of the driving vanishes, the multi- 
scale geometry of the critical percolation is dominant in providing the major trans- 
port characteristics for the DPRW lattice. The situation changes drastically when 
the strength of the driving increases above some level. With the systems departure 
away from the percolation point, the multi-scale features will soon be lost substi- 
tuted by the bulk-average nonlinearities. The fact that Eqs. ( [65] l and (|66]l above are 
formulated in terms of the system-average parameters, p and /, merely reflects that 
the system is allowed to appreciably depart from the state of marginal stability, or 
the SOC state (that means that pmax can be rather closer to 1 than to pc), and that 
the effect of overdriving readily calls for the global features to come into play. It 
is noted that, in general, the multi-scale properties due to SOC can coexist along 
with the global or coherent features; one example of this is substorm behavior of the 
dynamic magnetosphere 1451 . 

The end result of the discussion above is that the strength of the driving plays 
a crucial role in dictating both linear and nonlinear behaviors in the DPRW model. 
To obtain a pure SOC state the driving rate should go to zero fast enough as the 
critical point is approached. The main effect overdriving has on the DPRW dynam- 
ics is to excite unstable modes associated with periodic bursts in the particle loss 
current. Accordingly, the system auto-oscillates between a subcritical (pmin < Pc) 
and a supercritical (pmax > Pc) states in response to external forcing. The transition 
to auto-oscillatory dynamics signifies the increased role of global and nonlinear be- 
haviors in the strongly driven DPRW system as compared to a pure SOC system. 
The borderline between the two regimes corresponds to ~ {Pc — Pnim) The 
stability condition <C {pc — Pmm) /^+ has serious implications for the achievable 
SOC regimes. It poses one important restriction on the net deposition rate of the free 
particles against the longest relaxation time on the incipient percolation cluster. 



"Fishbone "-like instability 

To help judge the result obtained, let the critical exponents take their mean-field val- 



ues: z = 2, V = 1/2. In this limit Eqs. (65 i and (66 1 above reproduce, up to change 
of variables, Eqs. (13) and (14) of Ref. 1 134|. The latter set of equations appear in a 



Remember that, in the proposed model, the polarization charges act as conducting states for the 
motion of current-carrying particles. 
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basic theory of Alfven instabilities as a simplified model for the coupled kink-mode 
and trapped-particle system in a magnetically confined toroidal plasma where beams 
of energetic particles are injected at high power. The mode dubbed "fishbone" is 
characterized by large-amplitude, periodic bursts of magnetohydrodynamic (MHD) 
fluctuations, which are found to correlate with significant losses of energetic beam 
ions 1 1 34 1 . By comparing our Eqs. ( [65] l and ( [66| with Ref. II 13411 one can see that 
the lattice occupancy per site p corresponds to the effective resonant beam-particle 
normalized pressure within the q = I surface (here, q is the familiar safety factor 
used in tokamak research); pc corresponds to the mode excitation threshold; and 
the particle loss current / corresponds to the amplitude of fishbone. This direct cor- 
respondence between the two models suggests consider the instabiUty in Eqs. ( [65] l 
and (|66]) as analog "fishbone" instability for SOC dynamics. 

This correspondence is not really surprising. Mathematically, it stems from the 
resonant character of the fishbone excitation, implying that the energetic particle 
scattering process is directly proportional to the amplitude of fishbone 11331 11341 . 
This resonant property dictates a specific nonlinear twist to the fishbone cycle, dif- 
ferentiating it from other bursting instabilities in magnetically confined plasmas. It 
is this "resonant" twist observed in the DPRW model system that identifies the ana- 
log "fishbone" mode for SOC. We note in passing that the existence of an instability 
on the top of SOC dynamics conforms with the results of Ref. II135I in which the 
traditional (sandpile) SOC model has been modified by adding diffusivity, giving 
rise to periodic relaxation-type events as a function of the system drive, while a 
pure SOC state requires a vanishing drive 



The threshold character of fishbone excitation 

The fishbone belongs to a specific class of instabilities, the energetic particle modes 
(EPM), which appear in a magnetic confinement system when the energetic par- 
ticle pressure is comparable with the pressure of the thermal plasma 11331 11361 . 
The EPM constitutes a separate branch with a distinctive dispersion relation and its 
frequency and the growth rate crucially depend on the parameters of the energetic 
particle orbital motion. When the drive is strong enough, the EPM may be unstable 
despite having a frequency in the Alfven continuum where normal modes of the 
background plasma are typically strongly damped (the associated damping rate is 
proportional to the gradient of the phase velocity) 13711 . Instabilities in the Alfven 
continuum are often observed during intense neutral-beam injection II134|[T38 1, but 
they can also be excited by the energetic electrons generated experimentally by dif- 
ferent means: electron cyclotron resonance heating (ECRH) as on DIII-D tokamak 
II 1391 and lower hybrid (LH) power injection as in Frascati Tokamak Upgrade (FTU) 
experiments II 1401 . The phenomenon is thresholded in that it requires a critical level 
of the absorbed power The existence of the critical power, which we associate with 
the critical "driving" rate, Z+max, is well established experimentally (Fig. 8). 

Beside this share, the mode referred to in Ref. f 1351 is edge localized, non-resonant, fluid mode, 
in agreement with the diffusive nature of the added flux, but at contrast with the resonant mecha- 
nism of the fishbone excitation. 




42 



Alexander V. Milovanov 



> 



> 



1.5 
1.0 
0.5 
0.5 



-0.5 
11 
10 
9 



_ Plh llti^tm^Y^ 


1 ^fiilililiiJi 


1^' branch J 2"^ 


branch 


1 1 ECE Ch 9 


S, AHA.^ ECE Central 



0.20 



0.25 



0.30 
Time (s) 



0.35 



0.40 



Fig. 8 Plots of LH coupled power, fast electron temperature fluctuations, and central radiation 
temperature in FTU shot # 20865. During high power LH injection, an evident transition in the 
electron fishbone signature takes place from almost steady state nonlinear oscillations (fixed point; 
marked as P' branch) to regular bursting behavior (limit cycle; marked as 2"'' branch). The tran- 
sition is at f^H ~ 1-69 MW. It is noticed that the bursting behavior phase closely resembles that 
of well-known ion fishbones |134||138il and ECRH driven electron fishbones on DIII-D U39il . 
Adapted from Ref fT40l . 



The EPM dispersion relation i l31lll33lll34l ll361 when account is taken for the 
well-known "fluid" (this includes the background MHD and the energetic particle 
adiabatic and convective responses) and "kinetic" contributions due to the energetic- 
particle "compressions" can be obtained via asymptotic matching procedure, lead- 
ing, upon the fast and slow time scales are separated, to the frequency-dependent 
complex nonlinear parabolic equation, Eq. (15) of Ref. II131I . A remarkable fea- 
ture of this equation is that the nonlinearity due to the wave field is twisted with 
the free energy source term (here thought as the "driving" term). The main effect 
this twist has on the dynamics is that the EPMs are released in radially amplify- 
ing "avalanches" II133I (to visualize, think to a large mass of mud or snow moving 
rapidly downhill). This avalanching behavior which we associate with behavior of 
a strongly overly driven system in the presence of intense energetic particle popula- 
tion should be distinguished from the above "chain reactions of hopping motions," 
which are the avalanches in the DPRW SOC system at criticality. In the local limit, 
characterized by a Gaussian free energy source profile, Eq. (15) of Ref. 11311 is 
reduced to a complex NLSE, which is different from DANSE (46 1 in that it is dom- 
inated by the nonlinear properties, rather than by a competition with randomness. 
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Fractional nonlinear Schrddinger equation 

If one wants to go beyond the local limit, one may use a stretched Gaussian free 
energy profile instead II141II . In the latter case, the resulting equation is found to be 
a variant of fractional nonlinear Schrodinger equation, or FNLSE (presented here 
without derivation, see Ref. 1 13|) 

'^^^-^.,VLV?f (f,x) + C|f =0, (69) 

where 

V^f (r,x) - [\lr{\-q)]x-^^W{t,x) = [l/r(l -^)] J^yy\x~y\-"^{t,y) 

(70) 

is the so-called RieszAVeyl fractional derivative BtI ItTI . which is mathematically 
different from the Riemann-Liouville derivative in Eq. ( [35] l in that the integra- 
tion is performed through infinite limits; the operator VL^i is a generalization 
of the Laplacian; x denotes respective spatial coordinate; <7 is a fractional exponent 
(0 < <7 < 1) which corresponds with the exponent of the stretched Gaussian free en- 
ergy source profile |141 1; and is a normalization constant which carries the di- 
em^'' • s^' (we assume that h= 1). For attractive nonlinearity (here 



corresponding to ^ > 0), FNLSE ( 69 1 describes phenomena of self-delocalization 
of fractons (vibrational excitations of fractal networks) as well as associate phe- 
nomena of nonlocal soUtary waves ifTSll . Setting f (f,x) = \j/{x)exp{—i(Ot), from 



FNLSE (69 1 one arrives at the fractional envelope equation 

- i^,Vl,V^v^(x) + «v/(x) + CIrWI VW = 0, (71) 



which does not contain time differentiation. Remark that FNLSE ( |69| l is built on 
fractional derivatives in space, while time differentiation is conventional (integer). 



likewise to DANSE (46 1. This distinction between the fractional derivatives in space 
and time reflects the very special role 1131 11421 time plays in the dynamical equa- 
tions originating from quantum mechanics such as NLSE, by contrast with the ki- 
netic equations for transport and relaxation processes discussed above Ij47l l74l l75l . 



Mixed SOC-coherent behavior 

The idea of "fishbone" instability in self-organized critical dynamics is very ap- 
pealing as it addresses a type of behavior in which the multi-scale features due 
to SOC can coexist along with the global or coherent features. One example of 
this coexistence can be found in solar wind— magnetosphere interaction. Indeed it 
has been discussed by a few authors 11431 11441 11451 11461 that the coupled solar 
wind— magnetosphere— ionosphere system operates as an avalanching system and 
that there is a significant SOC component in the dynamics of magnetospheric storms 
and substorms, along with a coherent component 1 143 , 147 1 that evolves predictably 
through a sequence of clearly recognizable phases lil48J . Here, we advocate a way 
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of thinking lfT3l l99l in which the magnetospheric SOC component is associated 
with the properties of self-organization of electric currents and magnetic field fluc- 
tuations in the plasma sheet of the Earth's magnetotail (i.e., the "Sakura" model); 
whereas the coherent component is attributed to global instability of the cross-tail 
SOC current system and the phenomena of tail current disruption. The implication is 
that the dynamic magnetosphere survives through a mixed SOC-coherent behavior. 
In this spirit, we expect the input power due to magnetic reconnection at the Earth's 
dayside magnetopause to self-consistently control the departure from the state of 
marginal stability, or the SOC state, with stronger departures favoring the coherent 
features. By analogy with fishbone instability in magnetic confinement systems we 
suggest that behavior is thresholded in that there is a critical input power (critical 
reconnection rate at the dayside magnetopause), destabilizing the SOC component 
in the magnetotail. At this point, a portion of the cross-tail electric current will be 
redirected to ionosphere — thought as analog "ground circuit" (Fig. 9) — leading to 
a decrease in the tail current intensity (tail current disruption), and to a magneto- 
spheric disturbance, or a substorm. A model for the substorm cycle is obtained from 
the above coupled system of equations, Eq. ( [65| ) and ( [66| , where one identifies the 
driving rate, Z+, with the magnetic reconnection rate at the Earth's dayside mag- 
netopause; the lattice occupancy parameter, p, with the average normalized energy 
density of magnetic fluctuation field in the magnetotail current sheet; and the parti- 
cle loss current, /, with electric current in the ionosphere. We should stress that we 
consider a substorm as instability in the cross-tail SOC current system, which oc- 
curs on the top of self-organization to a critical state in the magnetotail, thus posing 
a coherent component over the dynamics. 



9 Phase transitions in SOC systems 

It was argued that the phenomena of magnetospheric substorm bear signatures en- 
abling to associate them with a second-order phase transition in the coupled per- 
turbation electric current and magnetic fluctuation system 1981 l99l . Indeed, when a 
current filament pops up in the plasma sheet, it interacts with the Harris-distributecT^ 



magnetic field in the lobes of the magnetotail (Fig. 9). The forces are such as to make 
the filament spontaneously change its orientation, and in fact a local minimum in the 
free energy profile occurs, which favors current disruption to reinstall stability. 

Let us address the phenomena of magnetospheric substorm from a more funda- 
mental perspective, namely, as part of the general problem of phase transitions in 
SOC systems. The main idea here is that some systems may spontaneously turn into 
a coherent state before they become SOC, since their evolution by itself drives these 
system to a competition between the SOC and coherent properties as a consequence 
of some nonlinear twist between associate order parameters. Other than substorms, 
this general approach may include phenomena like the L-H transition in magnetic 



According to Harris |1491, tlie dependence of the lobe field is given by hyperbolic tangent of the 
distance to the neutral plane. 




Percolation Models of Self-Organized Critical Phenomena 



45 




Fig. 9 Substorm in the Earth's magnetosphere. Red: The complex electric current system in the 
magnetotail current sheet. Blue: The magnetotail lobe field, B^^. A current filament popping up in 
the plasma sheet interacts with the Harris-distributed magnetic field in the lobes of the magneto- 
tail. The forces are such as to make the filament spontaneously change its orientation, favoring 
redirection of electric current to the ionosphere. 

confinement devices [1501, and the tokamaks as particular case, where the L-phase 
is associated with SOC |85|, and the H-phase is associated with spontaneously oc- 
curring coherent state. 

Subordination to SOC 

Let us consider a spatially extended system with some order parameter T, where the 
processes of self-organization develop a singularity at some value for f +°o. 
We assume that this singularity does not explicitly appear in the dynamics, implying 
that the system is, in this limit, critical and self-organized. The phenomena we are 
looking at appear when the system possesses a competing order parameter, which 
we shall denote by and for which Y acts as input control parameter The im- 
plication is that the dynamics of y/ is subordinated to the dynamics of Y via some 
intrinsic coupling mechanism. For simplicity, we shall assume that the order param- 
eter \j/ corresponds to a coherent behavior, which we envisage as competing with 
the emerging multi-scale features due to SOC. Thus, while the system is developing 
its singular (SOC) points, it may find it thermodynamically profitable to sponta- 
neously turn into the competing, coherent phase. As we shall see, this idea leads 
to a fractional extension of the Ginzburg-Landau equation |fT3]|151| , in which the 
conventional Laplacian is replaced by fractional Riesz/Weyl derivatives. 
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Generalized free energy expansion 

For the sake of mathematical convenience, we shall assume that the system ap- 
proaches its SOC point so closely, that the geometry is, to a good approximation, 
self-similar (fractal). In this regime, the distribution of the competing order parame- 
ter will be characterized by the diverging correlation length, owing to its coupling 
with r, so that both and T distributions are heavy-tailed. In a sense, the fractal 
distribution of T acts a fractal support for \//. The fact that the variation involves 
correlations on many spatial scales must have implications for the generalized form 
of the free energy expansion near the phase transition point, and in particular for 
the gradient term, where the usual assumptions of locality, permitting to write this 
term as a simple |VvV/(x)p, should be relaxed. Then a consistent generalization ac- 
counting for the integral effect of the correlations is obtained in terms of a Fourier 
convolution, V%\i/{x) = — q)]x^'^ -k\if{x), or the Riesz/Weyl fractional deriva- 

tive BtI ITTI . where q (0 < q < 1) characterizes the strength of spatial correlation 
and the local limit is reinstalled for ^ — > 1. It is this convolution which we expect to 
replace V^y/(x) when the gradient term is considered. Indeed, the following gener- 
alized free energy expansion in vicinity of the transition point holds lfT3lll51ll 



F = Fo+ dx 



4|V.Xx)|2+fl,|v/(x)|2 + lfo,|V^(x)r 



(72) 



where we have introduced three phenomenological expansion parameters a^, 
and bq, which may depend on the exponent q in general. 

Fractional Ginzburg -Landau equation 



Varying the functional in Eq. ( 72 1 over the complex conjugate (x) and consider- 
ing yf{x) and V^*(x) as independent order parameters, one obtains 

8F^ p J;c[^^V?V^(x)V?5v/*W+fl<,V^(-«)5v/*(x)+^<,|V/'(.x)|VW5v^*W] . 

(73) 

With use of the integration-by-parts formula IItTI 

<ix^l(x)V?^(x) = y dx(p2(x)V'^_^(p\{x) (74) 



equation ( 73 1 implies that 



5F= Jx[i4V^^VfvAW + fl^v^W + fe,|v/W|VW]5v^*W, (75) 



yielding, in view of the extremum 5F = 0, 

^^V^,V?v/W+fl^r(x)+fe,|V/(x)|VW =0. (76) 



Note that varying the integral in Eq. (72 1 over \//(x) leads to the conjugate equation 
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i^V^ ,V>*(x) +fl,V^*(x) +^7,1 v/WI V* W = 0, (77) 

which is physically identical to Eq. ([76]). Equation ( [76| has the mathematical struc- 
ture of the well-known Ginzburg-Landau equation il52l . in which the conven- 



tional Laplacian, V^, is generalized to V'L^V^. We consider Eq. (76i as fractional 
Ginzburg-Landau equation, or FGLE. One sees that FGLE appears as a natural tool 
in describing phase transitions in SOC systems, in much the same way as the con- 
ventional Ginzburg-Landau equation describes type II phase transitions in simpler 



systems. Remark that FGLE (76 1 is different from the fractional envelope equation. 



Eq. (71 1, in that it contains the opposite sign in front of Vl^Vf . 



The q exponent 

Likewise to traditional type II phase transitions one may argue that Uq changes 
sign at the critical point and that it linearly depends on variation of the input 
control parameter 11521 11531 . Then the subordination condition will imply that 
Qg = OLqiX ~ ^ ) for ^ ~^ with aq a constant which does not depend on T . 
Given that the system is driven so slowly that it develops a singular point as a 



result of self-organization, one predicts, with the use of FGLE ( 76 1, that the dis- 
tribution of the order parameter i// will be self-similar to comply with the scaling 
I ^'{x)\'^ o= x^^^, from which the Hausdorff dimension df — d--2q can be deduced. 
Focusing on the df value, because the order parameter is intrinsically coupled to 
r, one may expect that the and T distributions will be essentially the same and 
hence, when account is taken for the percolation fractal geometry of SOC, will be 
both characterized by the "hyperuniversal" relation, df — d ~ P /v (see section 1) 
where j3 and v are percolation critical exponents. The latter expression will be con- 



sistent with FGLE (76i, when q = j3/2v. This is the desired result. It shows that 
there is nontrivial fractional index of integro-differentiation in associate fractional 
Ginzburg-Landau equation and that behavior is nonlocal in general. Using known 
estimates for the parameters j3 and v [6, 7|, it is found that q = 5/96 for d = 2; 
q w 0.26 for d = 3; and q « 0.4 for d = 4. The mean-field value, holding for d >6, 
is <7 = L One sees that the mean-field case is local, as it should. Allowing for <7 — > 1 



in FGLE ( 76 1, one verifies that the conventional ("integer") Ginzburg-Landau equa- 
tion is readily reinstalled. By contrast, for ^ < 1, the dynamics are governed by 
an interplay between nonlocality and nonlinearity, likewise to the EPM excitation 



case, Eq. ( 69 1, and their correct description requires fractional extensions of respec- 
tively the Ginzburg-Landau and nonlinear Schrodinger equations consistently with 
the implication of the RieszAVeyl fractional operator. 



10 Overall summary and final remarks 



The concept of self-organized criticality, or SOC, proves to be a complementary tool 
in drawing a physical picture of the processes underlying the dynamics of systems 
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with many coupled degrees of freedom (i.e., the "complex" systems). In this report 
we have demonstrated the diversity of appropriate mathematical methods of describ- 
ing such processes, including fractional equations of the diffusion, relaxation, and 
Ginzburg-Landau (nonlinear Schrodinger) type, generalizing their standard counter- 
parts, as well as the formalism of discrete Anderson nonlinear Schrodinger equation 
(DANSE) — extending far beyond the usual scaling theories. Some connections to 
Hamiltonian models, paving the way to first-principle models of SOC phenomena, 
have been also discussed. These issues make the mathematical formalism of SOC 
an exciting and challenging problem. 

The main emphasis in the present work has been laid on percolation, recognized 
as a convenient and powerful framework in describing critical phenomena in com- 
plex systems. The percolation problem finds its significance in its relation with the 
fundamental topology (in terms of connectedness issues) and theory of fractional 
manifolds llT3l l23l . In this respect, some elements of conformal maps of fractals, 
along with the percolation problem on the Riemann sphere, have been addressed. 

To deal with dynamical problems involving feedback between the various de- 
grees of freedom such as the SOC problem we used the idea of random walks on a 
self-adjusting percolation cluster This idea is advantageous, as it makes it possible 
to employ the random walks in place of the usual lattice automated redistribution 
rules, and in this fashion to design a family of SOC models, that are more friendly 
from the standpoint of their analytical treatment as compared to their CA relatives. 
In fact, the random walks on percolation systems offer suitable analytical forms 
for the diffusion, charge-conduction, and the dynamic susceptibility properties, and 
their significance in the study of SOC phenomena can hardly be exaggerated. 

More so, we proposed a simple lattice model of self-organized criticality, the 
DPRW model, which addresses the SOC problem as a transport problem for electric 
charges (free particles and holes) on a dynamical geometry of the threshold perco- 
lation. The novel concepts of this model are: (i) a theory of self-organized criticality 
based on the analogy with dielectric-relaxation phenomena in self-adjusting random 
media, and (ii) prediction of a "resonant" instability of SOC due to the nonlinearities 
present. The system adjusts itself to remain at the critical point via the mechanisms 
of hole hopping associated with the random walk-like motion of lattice defects on a 
self-consistently evolving percolation cluster. 

With the random walk guide to lattice dynamics we could derive fractional 
analogs of the diffusion and relaxation equations, demonstrating the existence of 
multi-scale relaxation processes and of a broad distribution of durations of relax- 
ation events. In particular, we have shown that the relaxation to SOC of a slightly 
supercritical state is described by the Mittag-Leffler relaxation function, Eq. p2| ), 
and not by a simple exponential function as for standard relaxation. The ideal SOC 
state requires that the driving rate goes to zero faster than a certain scaling law as the 
percolation point is approached. The model belongs to the same universality class 
as the BTW sandpile, and should be distinguished from the DP-like SOC models. 

Thinking of holes as "excitations" of the marginally stable state, we considered a 
transport problem for the hole wave function in the context of DANSE equation with 
random potential on a lattice. An important feature which arises in this approach 
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is competiotion between nonlinearity and randomness. It was argued that above a 
certain critical strength of nonHnearity the Anderson locaHzation of the hole wave 
function was destroyed and unlimited subdiffusive spreading of the wave field along 
the lattice occurred. This subdiffusion is asymptotic |119|. We have seen that this 
problem of the critical spreading was intimately related with the outstanding prob- 
lem of transport along separatrices in large systems fl24]. With the recognition that 
the transition to unlimited spreading could be described as a percolation transition 
on a Cayley tree, a "self-organized" formulation of the phenomena of localization- 
delocalization in the presence of nonlinearity has been proposed. The results of this 
investigation have demonstrated the versatility of the DANSE formalism, which we 
believe to capture the essential key elements of self-organized critical phenomena, 
thus offering a general analytical framework for SOC. 

Overdriving the DPRW system near self-organized criticality was shown to have 
a destabilizing effect on the SOC state. The fundamental physics of this instability 
consists in the following. Because of rapid accumulation of the conducting sites, 
the system departs from the percolation point, and its geometry nonlinearly changes 
from fractal-like to crystalline-like. At this point, the conductivity of the system 
has greatly increased. As the lattice conducts more electricity, losses increase in the 
ground circuit. However, because the particle loss current has feedback on the lattice 
occupancy parameter, a cross-talk is excited between the systems average conduc- 
tivity response and the distance to the critical state. We have observed that the insta- 
bility cycle is qualitatively similar to the excitation of the internal kink ("fishbone") 
mode in tokamaks with high-power beam injection (the lattice occupancy per site p 
corresponds to the effective resonant beam-particle normalized pressure within the 
q= \ surface; pc corresponds to the mode excitation threshold; and the particle loss 
current / corresponds to the amplitude of fishbone). The instability is "resonant" 
in that the particle loss process is directly proportional to /. This resonant property 
dictates a specific nonlinear twist to the fishbone cycle, differentiating it from other 
bursting instabilities in magnetically confined plasmas. 

The excitation of "fishbone" instability in SOC systems leads to a type of behav- 
ior in which the multi-scale features due to SOC can coexist along with the global or 
coherent features (i.e., mixed SOC-coherent behavior). One example of this coexis- 
tence is found in the solar wind— magnetosphere interaction. We expect the concept 
of mixed SOC-coherent behavior be the plausible statistical picture for thresholded, 
dissipative, nonlinear dynamical systems in the parameter range of nonvanishing 
external forcing. In this respect, we suggest that some of the "extreme" events, or 
system-scale responses, observed in complex natural and social systems |154| may, 
in fact, be the fishbone-like instabilities of SOC predicted by the present theory. 

It has been shown that some systems may spontaneously turn into a coherent 
state before they become SOC, since their evolution by itself drives these system 
to a competition between SOC and coherent properties as a consequence of some 
nonlinear twist between associate order parameters. We discussed a generalized free 
energy expansion for a system with extended spatial degrees of freedom, in which 
the order parameter due to SOC acts as input control parameter for the competing 
coherent behavior Based on this expansion — which has involved the Riesz/Weyl 
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fractional operator in place of the standard (local) gradient — a fractional generaliza- 
tion of the well-known Ginzburg-Landau equation has been obtained variationally 
111511 . With the fractional Ginzburg-Landau equation, it should be possible to de- 
scribe the much observed L-H transition in magnetic confinement devices II150L as 
well as the phenomena of magnetospheric substorm and tail-current disruption. 

We believe it will be worthwhile to pursue the above considerations not only 
because they arise naturally in a basic theory we are considering, but also because 
questions of this kind might have feedback on seemingly very diverse phenom- 
ena beyond the specialized physics context of this study. We illustrate this on two 
examples from respectively finance and climate dynamics ll34l 11541 . Indeed col- 
lective phenomena in finance include as partial cases the bursting of speculative 
bubbles, market crashes, debt contamination (now-deepening in the euro-zone), and 
the spreading of bankruptcy and insolvency. A simplified toy-model here might be 
constructed as a variant of the DPRW model discussed above, with clusters of po- 
larization charges thought as asset market; holes thought as insolvency; and the 
phenomena of hole-hopping thought as debt spreading. In this context, we might 
predict that, when the market is de-regulated (i.e., the dynamics are random walk- 
like), periodic financial crises are virtually unavoidable and that the period between 
crashes is inversely proportional to the rate at which speculative (not absorbed by 
the real economy) capital flows into the market, thus "overheating" the system. 

Likewise to finance, many exciting questions such as the above arise in climate 
research and it would be of interest to study them from the more general perspective. 
The focus here is on existing periodic oscillations in Earth's global climate system, 
as for instance the El Nino/La Nina-Southern Oscillation, or ENSO, which is a two 
to seven year quasi-periodic climate pattern associated with warming (El Nino) and 
cooling (La Nina) of the ocean surface layer across the tropical eastern Pacific. The 
extremes of this oscillation — which is identifiable in the climate reconstructions 
since thousands of years — are blamed for the severe weather conditions affecting 
climate, habitats, and the economies in many regions of the world |155 |. In addi- 
tion, ENSO involves 1 156 1 interactions extending through different time scales with 
various climate phenomena, such as the seasonal cycle, interseasonal oscillations, 
and/or decadal oscillations 1 156111571 . 

It was argued 1 158 1 that ENSO is a naturally occurring, oscillatory mode of the 
coupled ocean-atmosphere system. We suggest that this mode is strongly driven in 
that it is excited when atmospheric forcing exceeds a certain critical level. Then 
the oscillation is of "fishbone" type, since the rules by which the mode is sustained 
will include features of unstable SOC dynamics |45|. One might expect that a sim- 
plified yet relevant toy-model here will take the form of two cross-talking vari- 
ables, Eqs. (65 I and (66 1, where p stands for air surface pressure and / stands for 
ocean surface temperature and that the period of the oscillation will be inversely 
proportional with atmospheric forcing strength. More advanced models might refer 
to FNLSE (69 1 with an account for cross-scale couplings between the various wave 
processes involved. Analysis of this general area remains to be carried out. 

In the same spirit, the phenomenon of fishbone might shed new light on glacial- 
interglacial climate changes, in particular, the ^ 100,000-year climate cycle, which 
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Thousands of Years Ago 

Fig. 10 Climate reconstruction over the past 420,000 years from the Vostok ice core, Antarc- 
tica. Blue: Temperature variation in degrees Celsius. Green: Carbon dioxide (CO2) concentra- 
tion in parts-per-million-by-volume, or p.p.m.v. Red: Dust concentration (dust concentrations are 
expressed in parts-per-million, or p.p.m., assuming that Antarctic dust has a density of 2,500 
kg-m^^). One sees that atmospheric concentrations of carbon dioxide correlate well with Antarctic 
air-temperature throughout the record, while anit-correlate with the dust concentration. Image and 
data credit: Ref fT59l . 

could be a globally induced unstable climate mode stemming from a cross-talk be- 
tween air-temperature and dust concentration (see Fig. 10). Support for this sugges- 
tion can be found in Antarctic ice records as discussed in Ref. II159I . The implication 
is that glaciations occur naturally through functioning of Earth's climate as complex 
system, thus being "inherently there" as the many degrees of freedom twist and 
couple with the exterior. All in all, the phenomenon of fishbone warns of repeating 
severe events being virtually unavoidable in driven systems. 

Generally speaking, the study of self-organized criticality phenomena is cur- 
rently transitioning from an emphasis on scaling and linear-response theories to an 
emphasis on understanding and predicting the nonlinear dynamics of systems with 
many coupled degrees of freedom. In many ways these tendencies are manifest in 
the introduction of DANSE, FNLSE, and FGLE equations discussed above. As this 
transition occurs, investigations — as much experimental as theoretical and numeri- 
cal — that pinpoint the nonlinear interactions in complex systems will increase in im- 
portance. In the geo-space plasma research, with the observations becoming multi- 
spacecraft and/or multi-point in scope, theoretical models are likewise to confront 
issues of nonlocality, self-organization, and build-up of correlations in the presence 
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Fig. 11 A cartoon illustrating the many co-going plasma processes in a toroidal magnetic con- 
finement system, viewed through the prism of complexity and self-organized criticality. A remake 
of Figure 1 from Ref. 1131 . In the bottom-right comer is an artist's view of the FAST tokamak. 
Courtesy of dr. Aldo Pizzuto. In the top-right comer is a "sandpile" watch: Time is running for 
fusion. The four satellites in the bottom centre represent the ROY project 11601 . 

of many co-existing plasma processes |fT3lll60ll . Similarly to geo-space exploration, 
research activities in fusion plasma are now arriving at a crucial juncture that neces- 
sitates the understanding of "complexity" in the accessible and relevant operation 
regimes of burning plasma. Indeed it is becoming clear that the important questions 
that will be receiving attention in the coming years, particularly with the devel- 
opment of ITER and DEMO scenarios, are addressed toward the comprehension 
of burning plasma state as being self-organized, thresholded, nonlinear dynamical 
system with many interacting degrees of freedom f l31l 11331 11411 L1611 . The ITER 
project is a major challenge on the way to controlled fusion burn. The specialized is- 
sues of complexity, nonlinear interactions, and SOC have found their significance in 
the recently formulated Fusion Advanced Studies Torus, or FAST, proposal II162L 
promoting an European Union satellite for ITER. An emergent way of thinking 
here is to recognize FAST as having a parallel in the geo-space exploration, the 
ROY mission concept [160] — a project in space research for a constellation of 
small, probe-like satellitesP£| aiming to investigate the dynamic magnetosphere as 



ROY means "Swarm" in Russian 
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complex system (see Fig. 11). We extrapolate that the cross-disciplinary effort of 
bringing these exciting projects to realize will open new avenues in the study of 
what proves to be one of the greatest theoretical challenges in the modem nonlinear 
physics, the paradigm of self-organized criticaUty. 



Table 2 Abbreviations used 


Abbreviation 


Expansion 


AE 


Alfv6n eigenmode 


AMPTE 


Active Magnetospheric Particle Tracer Explorers (satellite) 


AO 


Alexander-Orbach (conjecture) 


BTW 


Bak-Tang-Wiesenfeld 


CA 


Cellular Automation 


CTRW 


Continuous Time Random Walks 


ECRH 


Electron Cyclotron Resonance Heating 


ENSO 


El Niiio/La Nifia-Southem Oscillation 


EPM 


Energetic Particle Mode 


DANSE 


Discrete Anderson Nonlinear Schrodinger equation 


DEMO 


DEMOnstration Power Plant 


DIII-D 


DTTl-D fnknmak 


DP 


Directed Percolation 


DPRW 


Dynamic Polarization Random Walk 


ITER 


International Thermonuclear Experimental Reactor 


FAST 


Fusion Advanced Studies Torus 


FDE 


Fractional Diffusion equation 


FGLE 


Fractional Ginzburg-Landau equation 


FNLSE 


Fractional Nonlinear Schrodinger equation 


FTU 


Frascati Tokamak Upgrade 


KAM 


Kolmogorov-Arnold-Moser 


KWW 


Kohlrausch- Williams-Watts (relaxation function) 


LH 


lower hybrid (oscillation) 


L-H 


low-high (transition) 


MHD 


Magnetohydrodjmamic 


MW 


Megawatt 


NLSE 


Nonlinear Schrodinger equation 


ROY 


("Swarm") IKI-led multi-satellite geo-space mission concept 


SOC 


Self-Organized Criticality 
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